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Jnjtroduction 

In the last two decades, there have been extensive developments in computational 
aerodynamics, which constitutes a major part of the general area of computational fluid 
dynamics. Such developments are essential to advance the understanding of the physics of 
complex flows, to complement expensive wind-tunnel tests, and to reduce the overall 
design cost of an aircraft, particularly in the area of aeroelasticity. 

Aeroelasticity plays an important role in the design and development of aircraft, particularly 
modem aircraft, which tend to be more flexible. Several phenomena that can be dangerous 
and limit the performance of an aircraft occur because of the interaction of the flow with 
flexible components. For example, an aircraft with highly swept wings may experience 
vortex-induced aeroelastic oscillations. Also, undesirable aeroelastic phenomena due to the 
presence and movement of shock waves occur in the transonic range. Aeroelastically 
critical phenomena, such as a low transonic flutter speed, have been known to occur 
through limited wind-tunnel tests and flight tests. 

Aeroelastic tests require extensive cost and risk. An aeroelastic wind-tunnel experiment is 
an order of magnitude more expensive than a parallel experiment involving only 
aerodynamics. By complementing the wind-tunnel experiments with numerical simulations, 
the overall cost of the development of aircraft can be considerably reduced. In order to 
accurately compute aeroelastic phenomenon it is necessary to solve the unsteady 
Euler/Navier-Stokes equations simultaneously with the structural equations of motion. 
These equations accurately describe the flow phenomena for aeroelastic applications. 

At Ames a code, ENS AERO, is being developed for computing the unsteady aerodynamics 
and aeroelasticity of aircraft and it solves the Euler/Navier-Stokes equations. The purpose 
of this contract is to continue the algorithm enhancements of ENS AERO and to apply the 
code to complicated geometries. During the last year, the geometric capability of the code 
has been extended to simulate transonic flows, a wing with oscillating control surface. 
Single-grid and zonal approaches were tested. For the zonal approach, a new interpolation 
technique has been introduced. The key development of the algorithm was an interface 
treatment between moving zones for a control surface using the virtual-zone concept. This 
report summarizes briefly the work performed during the period, April 1, 1992 through 
March 31, 1993. Additional details on the various aspects of the study are given in the 
Appendices. 

V..-' 

Research Efforts 

The following specific objectives have been performed: 

1. The zonal procedure coded in ENS AERO was tested. Flows over a wing with 
oscillating control surface were simulated. The procedure was validated through 



comparisons with experimental data. The resultant code has been extended for wing-body 
configurations with oscillating control surface. See Appendix C and F. 

2. Parallel to the use of the existing zonal technique, more robust and generalized grid 
interpolation technique was investigated. A new procedure for zonal interface using the 
virtual zone technique has been added to ENSAERO. See Appendix D. 

3. Research for the turbulence model is undergoing. A problem with the Johnson-King 
model has been pointed out at Langley Research Center. Our independent calculation using 
their grid confirmed the observation. 

Results 

Virtual Zones 

Virtual zones are zones of zero thickness (for a finite volume formulation) which serve to 
transfer solid wall (or other) boundary conditions to an interface condition. Thus multiple 
boundary conditions can be imposed on a block face with the same flexibility as an 
interface condition. Virtual zones also decouple the process of volume grid zoning from the 
surface grid patches which define the aerodynamic configuration under study. Surface grid 
patches are required in order to impose the proper boundary conditions. The volume grid 
zones should be set up to obtain the proper mesh qualities required for numerical accuracy. 
Another advantage of the decoupling is that much fewer zones are now needed, thus easing 
the effort and time required to generate the grids about complex and realistic configurations. 

The orginal zoning capability of the ENSAERO code was extended by including the above 
capability of multiple interface conditions on a single block face. Since the code is a finite 
difference code the zones required an overlap at the boudaries of the zones to allow for the 
proper interblock communication (ie. interfacing). In this study a one cell overlap was 
chosen. For this kind of zoning the virtual zones of zero thickness used for the finite 
volume formulation is not appropriate. Instead the thickness of the virtual zones had to be 
expanded to include the extent of the overlap of the zones. In other words the virtual zones 
for the present formulation are now one cell thick. There is a slight mismatch of a half cell 
thickness between the location of the actual solid wall and the location where the virtual 
zones applies the solid wall boundary condition. This mismatch does not occur with the 
finite volume formulation. However in the present case the slight mismatch had no 
discernible influence on the overall flow field, especially in the case where flows at the 
ends of the flaps and wings are treated as viscous, rather than as inviscid. 

The test case considered was a clipped delta wing with an oscillating trailing edge control 
surface. The wing planform is shown in Fig. 1. The wing has a leading edge sweep angle 

of 50.4 deg and a 6 % thick circular arc airfoil section. At M„ = 0. 9 and a = 3 deg, both a 
leading edge vortex and a shock wave are present on the upper surface of the wing. The C- 
H grids of the three zones consist of 151 x 13 x 34, 151 x 15 x 34, and 151 x 20 x 34 
points from the inboard to the outboard zones. Since the experiment was conducted using a 
Freon test medium, the ratio of specific heats, y> was set to 1.135 in the present 
computations. 

Figure 2 shows the unsteady pressures with the control surface oscillating at a frequency of 

8 Hz and an amplitude of 6.65 deg at = 0.9 and a = 3 deg and Re c = 17 x 10 6 based 
on the root chord. The results are shown as the amplitude and phase angle of the upper 
surface at the three span stations as indicated in the figure. In general, the agreement with 
the experimental results is good. Because the accuracy of experiment had its own 
limitations, the virtual zone results are also compared with the single grid results. As 



shown there is quite a discrepancy between the virtual zone results and the single grid 
results (the 15 1 x 44 x 34 grid), especially in the amplitudes at the center of the control 
surface. The discrepancy is however most likely due to the gap that was introduced 
between the flap and wing in the single grid case to accommodate the shearing grid. If the 
gap is reduced by increasing the spanwise resolution of the wing and control surface, the 
computational results of the refined single grid case (151 x 87 x 34) approach those of the 
virtual zone. This particular example demonstrates the importance of simulating the 
geometry of the control surface/wing configuration accurately. 

Wing-Body Configuration with Oscillating Control Surfaces 

Next the code has been extended for unsteady Navier-Stokes simulations of transonic 
flows over a rigid arrow-wing body configuration of supersonic transport-type aircraft with 
oscillating control surfaces. The H-H topology grid was used. The ICEM DDN CAD 
software system by CDC was used to generate the surface grid. Then the volume grid was 
generated by using HYPGEN code. To treat the control surface movement, the single-grid 
approach was used. Figure 3 shows the geometry of the wind tunnel model. The 
configuration is a thin, low aspect ratio, highly swept wing mounted below the centerline 
of a slender body. The wing is flat with a rounded leading edge. The body is extended to 
downstream. This half-span grid consists of 110 points in the streamwise direction, 116 
points in the spanwise direction, and 40 points normal to the body surface, in total of 
510,400 points. In the following computations, the grid is further divided into the upper 
and lower grids at the wing and the H-topology cut condition is provided through a zonal 
interface. For the full-span configuration, the grid is mirrored to the other side and thus the 
number of the grid points is doubled. It should be noted that the exact wing tip definition 
was not available so the tip thickness was decreased to zero across three grid points. 

Figure 4 shows the steady pressures compared with experiment at three spanwise sections 

for the half-span configuration. The flow conditions consists of = 0.85, a = 8 deg, 8 = 

0 deg and Re s = 9.5 x 10 6 based on the mean aerodynamic chord. The leading-edge vortex 
is captured well by the computed results. Figure 5 shows the corresponding case with a 
flap deflection of 8 = 8 deg. The effect of the deflection is apparent at the outboard sections. 
Figure 6 illustrates the instantaneous antisymmetric position of oscillating control surfaces. 
Comparison of response characteristics between symmetric and antisymmetric control 
surface motions on the right and left wings is being studied. 

Concluding Remarks 

The geometric capability of the code, ENSAERO, has been extended to wing-body 
configuration with oscillating control surface. Comparisons with available experimental 
data show good agreement 

The future research plan is to extend the code toward a complete aircraft The next mile- 
stone is Navier-Stokes computation with engine thrust. For the code validation, good 
experimental data will be required. 

Research in turbulence modeling is ongoing. The models we tested often show dependency 
on the numerical dissipation. Further improvements will be studied for reliability. 




Fig. 1 Planform of clipped delta wing with tailing edge flap. 
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Fig. 2 Comparison of unsteady pressures between virtual-zone and single-grid 
computations with experiment. 
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Fig. 3 Wind-tunnel model geometry of an arrow-wing body configuration. 






















Fig. 5 Comparison of computed steady pressures with experiment. 






Fig. 6 Instantaneous flow field with antisymmetrically oscillating control si 
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Introduction 

B ODY-CONFORMING coordinate transformations of a 
fluid conservation law are generally used in computa- 
tional fluid dynamics. The associated metrics must satisfy 
certain geometric identities to maintain the global conserva- 
tion for numerical solutions. 1 These metrics are called free- 
stream capturing (or preserving) metrics. Numerical tech- 
niques are known to capture the freestream on stationary 
grids. 2 " 4 However, the extension of such a formulation for 
moving grids is not straightforward. The error introduced in 
forming the time metrics has been overlooked because it is 
negligible in most cases, but it can be significant in certain 
applications such as helicopter rotor flows. 5 Rigorous formu- 
lations based on the types of grid motions were discussed in 
Ref. 1, and demonstrated, for example, in Ref. 6. The present 
study describes detailed formulas that can be used in both 
finite volume (FV) and finite difference (FD) methods for 
constructing freestream capturing metrics in space and time. 


Finite Volume Formulation 

Geometric Identities and Freestream Capturing 

The integral form of a conservation law for a given cell can 
be written as 

[ Q dV-\ Q dK+'i d> n F dS dt =0 (1) 

J k(/ 2 ) J ku,) J /, J s(/> 

where V(t) is the cell volume and ndS(t) is a vector element of 
surface area with outwardly normal n . Considering the Euler 
equations, Q is a vector of conserved variables, viz., density, 
momentum, and energy, and F is the flux tensor of Q. The 
flux Fean be decomposed into the flux in the stationary frame 
F st and the contribution due to surface element velocity v as 


F = F n - vQ 


( 2 ) 
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For the freestream condition 0® and the geometric 
identities can be derived from Eq. (1) as follows: 


•/2 r* 

Qo>[ V(t 2 ) - V(ti)) =-F s t ® • O ndS dt 

J'i Js 


*/2 r 

+ Q® O rt • v dS dt 
J / 1 J s 


( 3 ) 


On a stationary grid, the first geometric identity can be ob- 
tained as 


To satisfy Eq. (5) numerically, it is essential to compute the 
right-hand side correctly. It can be rewritten as 

I (D n • v dS dt = 1 (j^Sn * v) dr = f Sn • v d/") (10) 
J/| Js J/.Vell / cell \ J,, / 

Let V s be a volume swept by a surface S between the time 
interval [f If / 2 ]: 

V s = I Sn * v dt (1 1) 

J'i 

Let S(r,) = S, 5 6 2 and S(t 2 ) - S r The volume V s can be 
computed similar to Eq. (8) as 

Ks 1263 = ^122'1'566'S' 

= ^(Sii's's + ^ 122 ' r + 5, 562 ) • (r 6 - - r,) (12) 

where S ir 5 ' 5 and 5, 22 ' r are the surface vectors in space and 
time domain. Then, Eq. (5) is satisfied numerically as 


v(h) - v{t t ) = '£v s (13) 

cell 

This formula requires the position of the grid at r 2 , not the grid 
velocity v itself as per Eq. (1 1). The position of the grid at the 
next time level must be given explicitly before computing the 
fluid motion. For aeroelastic simulations, for example, the 
structural equation of motion has to be coupled with the fluid 
equations only at explicit time levels. A fully coupled implicit 
formulation requires an iterative approach between the fluid 
and structural equations. 


(j) n dS = 0 (4) 

This gives a mathematical expression for a closed cell, which is 
a requirement for any grid system. The corresponding time 
integration in Eq. (3) will be satisfied for moving grids by 
applying Eq. (4) at each time level. From the rest of Eq. (3), 
the second geometric identity can be obtained as 

V(t 2 ) - V(t x ) = 0) n - v dS dt (5) 

J / 1 J S(/) 

This equation is essential for moving grids because it repre- 
sents the conservation of volume for a moving cell. To guaran- 
tee the global conservation for numerical solutions, these geo- 
metric identities must be satisfied numerically. 

Figure 1 shows a regular hexahedral cell. We assume that all 
edges are straight lines. Let r be the position vector of a point 
in space. The formulas for the surface vector S (note that S is 
taken in the positive coordinate direction here) and the volume 
V can be defined as 


Geometric Identities for Time-Differential Form 

The time-differential form of Eq. (1) is often used as a basis 
of the conservation law. Starting with the time-differential 
form, one can obtain similar geometric identities as in the 
previous section. However, this approach introduces an addi- 
tional geometric identity. 

The first geometric identity derived from the time-differen- 
tial form is identical to Eq. (4). Then, the rest of the equation 
can be written as 

dV I 

— = CD n ■ v dS (14) 

^ Js 

To proceed with the surface integration on the right-hand side, 
let us break down the surface element velocity v into its 
components. Let r and r 0 (/) be the position vectors of a point 
in space and the origin of the noninertial frame, respectively. 
Also let v 0 (O and Q(f) be the velocity and angular velocity of 
the noninertial frame relative to the inertial frame, respec- 
tively. These velocities represent the rigid motion of the grid. 
In general, 


^1562 = l / 2 (r 6 - r,) x (r 5 - r 2 ) (6) 

‘5*1562 = (^56 - **12) X (05 “ r 2t) (7) 

V\ 2345678 ~ ^ ( 5,485 + 5 , 2 3 4 + 5 i 562 ) * — T ,) ( 8 ) 

where r 56 = Vi(r s + r 6 ), and so on (see Ref. 1, for example). 
Note that Eqs. (6) and (7) result in the same expression. In 
fact, the surface vector is defined uniquely as long as the edges 
are straight lines. Also note that there are other consistent 
ways to compute the cell volume. However, Eq. (8) is the 
simplest form. 

Let 5,562 = S l562 /f. With either Eq. (6) or (7), Eq. (4) is 
satisfied numerically as 

ES/I = 0 (9) 

cell 


V = v 0 (O + fi(r) x [r - r 0 (t )] + v c (15) 

where v c (/) is the surface element velocity relative to the 
noninertial frame, which leads to change of a cell in time. 
From Eqs. (14) and (15), one obtains 

= (v 0 - fl x r 0 ) • ^<j) n dS^ + fi ■ ^<j) r x n dS^ 

+ ^<j) n • v c dS^ (16) 

The first integral on the right-hand side is zero due to the first 
geometric identity, Eq. (4). The second integral represents the 
rigid rotation of grid and does not contribute to the change of 
the cell volume. The remaining terms give the second geomet- 
ric identity similar to Eq. (5). 
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The third geometric identity can be obtained from the sec- 
ond integral in Eq. (16) as 


r x n dS = 0 


(17) 


aging procedure to numerically satisfy the differential chain 
rule. However, it is not straightforward to satisfy the dis- 
cretized form of the GCL. The discretized form of the GCL 
can be written as 

= M5 s c(S* * r T ) + 5,(5* • r T ) + 6 f (5' * r 7 )] (24) 


The discretized form of the third geometric identity can be 
expressed as 


£5rxn=0 (18) 

cell 

This is the additional requirement, which is concealed unless 
the grid movement is broken down properly. It must be satis- 
fied whenever the grid rotates. Reference 1 introduces the area 
moment 

M = J r x n dS 

to satisfy the discretized geometric identity, Eq. (18), on the 
hexahedron: 


^1562 = ri 65 x 5 i 6 5 + r 126 X *$126 (19) 

where r l65 = x A(r x + r 6 + r 5 ), 5 165 = Vi(r 6 - rO x (r 5 - r,), and 
so on. Note that M X562 * r 1562 x 5 i 562 . The expression r 1562 x 
5 i 562 is not well defined for computing area moment because it 
results in a nonzero sum over a cell. In contrast, Eq. (19) is 
well defined so that the sum over a cell is zero. The nonzero 
error introduced by the use of the inconsistent area moment 
has a unique feature. A simple analysis 7 shows that the error 
disappears on a hexahedron with parallelogram surfaces, for 
example, on the Cartesian grid. For a hexahedron with arbi- 
trary quadrilateral surfaces, the error remains constant rela- 
tive to the cell volume. This error may be ignored only if the 
effect of the Coriolis force is negligible. For example, the error 
is small for a flow over a wing oscillating at small amplitude. 
Note that the use of the surface area vector for evaluating the 
surface moment vector becomes consistent on a triangular 
surface. The use of the tetrahedral cell allows the most com- 
pact and consistent metric formulation, although it results in 
unstructured-grid formulations. 

Finite Difference Formulation 

The differential form of Eq. (1) has been used widely in FD 
formulations. With a generalized coordinate transformation, 
r = r($, 7 /,f,r) and t = r, the FD metric terms can be expressed 
as 


y=5«.f,U r -f,xr f = S« (20) 

and 


( 21 ) 

where r T = v and the transformation Jacobian 7= l/V. 
Analogous definitions can be derived for the other directions. 
The cell volume as defined by Eq, (8) is different from the 
inverse of the transformation Jacobian in a discretized form. 
Nevertheless, the cell volume can be applied to the FD method 
with a scaling factor of one-eighth. The differential forms of 
the geometric identities are 

(■$«)« + (S'), + (5 D r =0 (22) 

V r = (5* * r r \ + (5* * r r ), + (S> • r r ) f (23) 

Equation (23) is the differential statement of the geometric 
conservation law (GCL). 8 

The discretized form of Eq. (22) can be satisfied by the 
consistently differenced metrics, 3 which are based on an aver- 


where 6 indicates the central-difference operator in each coor- 
dinate direction. Analogous to the derivation of Eq. (18), if 
the grid is moved in a rigid rotation, V T = 0 and r T = Q x r. 
Then the left-hand side of Eq. (24) equals zero. However, the 
right-hand side becomes 5*(r x 5 s ) + 5,(r x S v ) + 5£r x 50 * 
0. Thus, the freestream will not be captured by solving Eq. 
(24) with the FD method. The GCL implies only a necessary 
condition to capture the freestream, not a sufficient condition. 

The analysis of the FD formulation can be simplified with 
the aid of the previous discussions for the FV formulations. 
Following Ref. 1, let the edges of the hexahedron in Fig. 1 be 
redefined as a double-sized cell in the FD grid with r x = 
r,_ ij_ \jc _ i, r 2 = r,*_ u+ i t *_ i, . . . , **3 = r i+ ij+ i,*- 1 » r s = 
r i+ u_ l%k + t. Also, let the time level advance from to t 2 . A 
simple calculation 7 shows that Eq. (7) is equivalent to the 
consistently differenced metrics. 3 Thus, the surface vector 
evaluations, Eqs. (6) and (7), on the double-sized cell can be 
regarded as the evaluations of the freestream capturing space 
metrics for the FD method. The FD time-metric evaluation is 
also obtained from the FV method. Time integration of Eq. 
(24) from r, to t 2 easily results in Eq. (13). Thus, for example 
in the £ direction, replacing S * * r T in the right-hand side of Eq. 
(21) with the time average, K s */A/ of Eq. (11), the freestream 
capturing time metrics can be obtained as 


1 / 

J 


\_ 

At 


S* • r T d t 


Vsi 

At 


(25) 


Note that £, defined here is costly but contains all of the 
information about the movement of a cell surface, such as 
translation, rotation, and deformation. In contrast, Eq. (21) is 
a simple product of surface area and velocity of cell centroid 
and can represent only a translational motion. 

The freestream subtraction technique will be useful for the 
rigid motion of the grid instead of the rigorous, costly evalua- 
tion of Eq. (25). Note that the subtraction is required only for 
the time-metric terms with the use of the freestream capturing 
metrics in space. 5 


Concluding Remarks 

This Note discusses the freestream capturing and the geo- 
metric identities for a fluid conservation law. To guarantee 
global conservation in numerical solutions, certain geometric 
identities must be satisfied. Based on the full integral form of 
the conservation law, Eqs. (7), (8), and (13) will guarantee the 
freestream capturing. Based on the time-differential form, 
another condition [Eq. (18)] must also be satisfied. However, 
when the differential form is used, the global conservation is 
not trivial. Considering an FV cell on the FD grid, the 
freestream capturing metrics in space and time can be con- 
structed from the FV formulations. Such an approach for 
evaluating the time metrics is costly but guarantees the global 
conservation for an arbitrary motion of the grid. 
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Unsteady Shock- Vortex Interaction on a Flexible Delta Wing 

Shigeru Obayashi* and Guru P. Guruswamyt 
NASA Ames Research Center , Moffett Field, California 94035 

Unsteady Navier-Stokes computations have been carried out for simulating transonic flows over a clipped 
delta wing undergoing oscillatory and ramp motions, including flexibility. The implicit upwind algorithm has 
been validated by comparing the solutions with experimental data for the oscillatory pitching motion cases. The 
numerical and experimental results agree well at moderate angles of attack, where a leading-edge vortex develops. 

The ramp motion cases have demonstrated the effects of unsteadiness of the flowfield and structural flexibility 
on the wing responses. For the 10-deg ramp motion, a vortex breakdown is observed. The inviscid interaction 
with the shock wave plays an essential role in the process of the breakdown observed in the present calculation. 


Introduction 

I N the last decade, there have been extensive developments 
in computational fluid dynamics toward the prediction of 
the three-dimensional flowfield about complex geometries at 
moderate and high angles of attack. Among the characteristics 
of flows over aircraft, the behavior of the flow over delta 
wings is of strong interest for high-speed aircraft because of 
the nonlinear lift increase due to the leading-edge vortex. 
Effects of flexibility can further influence the nature of flows 
on such wings. Steady-state flow problems associated with 
delta wings have been widely investigated computationally 
(for example, Refs. 1-3). Several advanced studies have also 
been performed for unsteady vortical flow calculations at sub- 
sonic and supersonic Mach numbers. For example, Ref. 4 
presented a conical flow computation on a rigid wing and Ref. 
5 presented a subsonic three-dimensional computation on a 
flexible wing, both involving vortical flows. 

Numerical methods can play an important role in comple- 
menting expensive wind-tunnel tests, particularly in the area 
of aeroelasticity. An aeroelastic wind-tunnel experiment is an 
order of magnitude more expensive than a similar rigid-body 
experiment involving only aerodynamics. By complementing 
the experiments with numerical simulations, the overall cost 
of the development of aircraft can be considerably reduced. 
Thus development of a numerical method is desired for sim- 
ulating aeroelastic phenomenon. To accurately compute 
aeroelastic phenomena it is necessary to solve the unsteady 
Euler/Navier-Stokes equations simultaneously with the struc- 
tural equations of motion. 

Recently a code, ENSAERO, was developed to compute 
aeroelastic responses by simultaneously integrating the Euler/ 
Navier-Stokes equations and the modal structural equations 
of motion using aeroelastically adaptive dynamic grids . 6 7 An 
upwind algorithm was implemented into the code and the 
resulting code was successfully applied to compute transonic 
flows over a typical fighter-type wing undergoing oscillatory 
motion . 8 
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The purpose of this article is to examine the capability of 
the present numerical method by simulating unsteady tran- 
sonic flows on a clipped delta wing, which contain a leading- 
edge separation. So far, unsteady transonic flow computations 
have dealt with the motion of shock waves on wings only at 
small angles of attack. The present study focuses on the for- 
mation and motion of the leading-edge vortex on a clipped 
delta wing at transonic speeds as well as motion of the shock 
wave. The code validation has been done for steady and un- 
steady pitching cases by comparing surface pressures with the 
experimental data . 9 Computations for ramp motion cases have 
demonstrated effects of the unsteadiness of the flowfield and 
the structural flexibility on the aeroelastic responses. 

Governing Aerodynamic Equations 
and Approximations 

The nondimensionalized, thin-layer Navier-Stokes equa- 
tions used in this study can be written in conservation-law 
form in a generalized body-conforming curvilinear coordinate 
system for three dimensions as follows: 

d r Q + d ( £ + d v F + ^G = — d ( G v (1) 

where r = t, £ = £(x,y,z,t), 17 = 17 (x,y,z,t), and f = l(x,y,z,t). 
In the present article, the £ and 17 directions are along the 
streamwise and spanwise directions of a wing, respectively. 
The viscous derivatives associated with these directions are 
dropped. In contrast, the f direction is normal to the wing 
surface, and thus the viscous derivatives are retained. 

The vector of conserved quantities Q and the inviscid flux 
vector F are 

-pi r pv ~ 

1 pw l put + T) x p 

Q = J pv , F = - pvV + T yp (2a) 

pw pwV + 7)zP 

e pHV - r\ t p 

where H is the total enthalpy and V is the contravariant ve- 
locity component. The time metric is related to the grid ve- 
locity as 

Vt = ~V x x t - y y y t - 17 Z z t ( 2 b) 

The Cartesian velocity components w, v, and w are nondi- 
mensionalized by the freestream speed of sound a the den- 
sity p is nondimensionalized by the freestream density p M ; the 
total energy per unit volume^ e is nondimensionalized by pjH. 
For the £ and f directions, E and G can be defined similarly. 
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= \ (“ 2 + + w\ + pr{y l _ 1} (a\ 

where Re is the Reynolds number, Pr is the Prandtl number, 
a is the speed of sound, and J is the transformation Jacobian. 
Pressure is related to the conservative flow variables Q through 
the equation of state for a perfect gas 

p = (r - 1) {e - \ (« 2 + v 2 + ^)} (3) 

where p is the fluid density and e is total energy per unit of 
volume of the fluid. See Ref. 7 for detailed definitions. 

The viscosity coefficient p in G v is computed as the sum of 
/x, + ix, where the laminar viscosity p, is taken from the 
freestream laminar viscosity, assumed to be constant for tran- 
sonic flows, and the turbulent viscosity p ; is evaluated by the 
Baldwin-Lomax algebraic eddy-viscosity model. 10 Since the 
flowfield to be considered in this article contains a leading- 
edge separation, it is important to apply the modification of 

Circular-arc airfoil 
t/c = 0.06 

L.E. sweep angle = 50.4° 

Area = 1635.88 in 2 
Span = 45.08 in. 

Root chord = 63.55 in. 

Tip chord = 9.03 in. 

Taper ratio = 0.1421 



Fig. 1 Planforir °«*ometry of clipped delta wing and typical flow 
structure. 


Fig. 2 Clipped delta wing and grid distributions at the root section. 



Fig. 3 Comparison of computed steady pressures using the standard 
and modified Baldwin-Lomax turbulence models with experiment, M „ 
= 0.90, a = 3.97, Re c = 17.6 x 10 6 . , Modified; , un- 

modified; □, A, experiment, NASA TP-2594. 

the turbulence model originally developed for crossflow sep- 
aration by Degani and Schiff. 11 The modification improves 
the pressure prediction as shown later, even though the sep- 
aration is fixed at the sharp leading edge in the following 
applications. 

Although the nature of interaction between vortex and shock 
wave is predominantly inviscid, the viscous terms are impor- 
tant to compute right vorticity. For example, a test calculation 
using the Euler equations showed that the leading-edge vortex 
is weaker because the inviscid model does not resolve the 
shear layer properly. 

Numerical Solution Procedure 

Among upwind algorithms, a streamwise upwind algorithm 
has recently been developed and applied to steady-state prob- 
lems of transonic flows over wings 12 and vortical flows over 
a delta wing 13 on fixed grids. Most multidimensional upwind 
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algorithms are first constructed in one dimension and then 
extended to multiple dimensions by applying the one-dimen- 
sional procedure in each coordinate direction. On the other 
hand, the present method uses the local stream direction, flow 
velocity, and pressure gradient to construct the upwinding. 
The switching of flux evaluations always takes place at sonic 
values, where transonic shock waves may be located. There- 
fore, this method follows the flow physics more closely than 




Fig. 4 Comparison of computed unsteady pressures using the stand- 
ard and modified Baldwin-Lomax turbulence models with experiment, 
Af x = 0.90, a m = 3.97, d = 0.46, Re c = 17.6 x 10 6 , k = 0.5919. 
, Modified; , unmodified; □, experiment, NASA TP-2594. 



a) Mode 1 : Frequency = 6 Hz 



Chord 

c) Mode 3: Frequency = 18 Hz 


Fig. 5 First four mode shapes 


the coordinate upwind methods. The computed results con- 
firmed the higher resolution of the present algorithm over the 
central-difference method as well as over other upwind meth- 
ods. 13 In this work, the streamwise upwind algorithm is ap- 
plied to compute the inviscid cell-interface fluxes. A second- 
order central-difference evaluation is applied to the viscous 
term. The complete algorithm can be found in Refs. 8 and 13. 

An implicit method is used for the time integration because 
the computational efficiency of the method is critical for ex- 
pensive unsteady viscous calculations. The method chosen 
here is the LU-ADI (lower-upper factored, alternating di- 
rection implicit) method 14 because it requires only scalar bi- 
diagonal matrix inversions. See Refs. 15 and 16 for additional 
details. 

Aeroelastic Equations of Motion 

The governing aeroelastic equations of motion of a flexible 
wing are solved using the Rayleigh-Ritz method. In this method, 
the resulting aeroelastic displacements at any time are ex- 
pressed as a function of a finite set of assumed modes. The 
contribution of each assumed mode to the total motion is 
derived by Lagrange’s equation. Futhermore, it is assumed 
that the deformation of the continuous wing structure can be 
represented by deflections at a set of discrete points. This 
assumption facilitates the use of discrete structural data, such 
as the modal vector, the modal stiffness matrix, and the modal 
mass matrix. These can be generated from a finite element 
analysis or from experimental influence-coefficient measure- 
ments. In this study, the finite element method is used to 
obtain the modal data. 

It is assumed that the deformed shape of the wing can be 
represented by a set of discrete displacements at selected 
nodes. From the modal analysis, the displacement vector {</} 
can be expressed as 

{d} = [<t>]{ q} ( 4 ) 

where [</>] is the modal matrix and {^} is the generalized dis- 
placement vector. The final matrix form of the aeroelastic 
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b) Mode 2: Frequency = 8 Hz 



Chord 

d) Mode 4: Frequency = 28 Hz 


frequencies of clipped delta wing. 
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Fig. 6 Comparison of sectional lift responses between rigid and flex- 
ible wings in 4-deg ramp motion, M * = 0.90, Re c = 15.0 x 10 6 , A 
= 0.04. , Flexible; , rigid; — • — , a (rigid + elastic). 


equations of motion is 

[A*M + [GM + [*]{,} = {F} (5) 

where [A/], [G], and [£] are modal mass, damping, and stiff- 
ness matrices, respectively. {F} is the aerodynamic force vec- 
tor defined as (typUll^^A^ACp} and [A] is the diagonal 
area matrix of the aerodynamic control points. 

The aeroelastic equation of motion [Eq. (5)] is solved by 
a numerical integration technique based on the linear accel- 
eration method. 17 

Aeroelastic Configuration Adaptive Grids 

One of the major deficiencies in computational aerody- 
namics using the Navier-Stokes equations lies in the area of 
grid generation. For steady flows, the advance techniques 
such as zonal grids 18 are being used. Grid generation tech- 
niques for aeroelastic calculations, which involve moving com- 
ponents, are in early stages of development. In Ref. 7, aero- 
elastic configuration adaptive dynamic grids were successfully 
used for computing time-accurate aeroelastic responses of swept 
wings. In this work, a similar technique is used. 

Results 

Numerical schemes used for flow calculations in aeroelas- 
ticity must guarantee the correct calculation of amplitude and 
phase of unsteady pressures. To verify the accuracy of the 
present code for simulating the complicated flowfield con- 
taining a leading-edge vortex and a shock wave, test cases are 
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Time (sec) 

Fig. 9 Comparsion of sectional lift responses between rigid and flex- 
ible wings in 10-deg ramp motion, M x = 0.90, Re c = 15.0 x 10 6 , A 
- 0.04. , Flexible; , rigid; — • — , a (rigid + elastic). 
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Chord 

Fig. 10 Deformation of flexible wing at t - 0.107 s in 10-deg ramp 
motion, = 0.90, Re c = 15.0 x 10 6 , A = 0.04. 


chosen from the experiment on a clipped delta wing undergo- 
ing prescribed pitching motion . 9 Since the experiment was 
conducted using a Freon test medium, the ratio of specific 
heats y is set to 1.135 in the following computations. 

Steady Pressures 

Steady-state calculations have been performed to examine 
the validity of the numerical procedure and the computational 
grid. The model planform geometry is shown in Fig. 1. The 
wing has a leading-edge sweep angle of 50.4 deg and a 6 %- 
thick circular-arc airfoil section. A typical flow structure is 
also illustrated in the figure. For this planform with a highly 
swept leading edge, a strong leading-edge vortex can form at 
moderate angles of attack. This vortex can interact with a 
shock wave at transonic Mach numbers. The lines. A, B, and 


C, indicate the spanwise locations of the pressure orifices in 
the experiment. Figure 2 shows the grid generated algebrai- 
cally in the C-H topology. The 17 , and £ coordinates rep- 
resent the chordwise, spanwise, and normal (to the wing sur- 
face) directions, respectively. The grid contains 151 x 25 x 
34 points. 

The present grid has less number of grid points when com- 
pared with the grids used in the typical steady-state compu- 
tations . 1-3 The number of grid points was determined to com- 
promise the accuracy and the total computational time. 
Unsteady computations require an order more computational 
time than the steady-state computations. With the present 
grid, a typical unsteady case can be computed within 5 h by 
using a Cray YMP computer with a single processor. (The 
code requires about 19 /is per grid point per time step for a 
flexible-wing case.) Then, three different time-step sizes were 
tried for each case to check the dependency on time-step sizes. 
Thus, in total, it took about 17 h to complete one unsteady 
pitching case including the corresponding steady-state com- 
putation. 

The present grid was chosen from several candidate grids 
by comparing the pressure distributions computed using each 
grid with the experiment for a typical flow condition contain- 
ing both a leading-edge separation and a shock wave at steady 
state. Among those grids, the H-O grid topology was also 
considered. However, the results showed that it did not have 
any advantage over the C-H grid topology either in capturing 
a leading-edge vortex (because of the moderate sweep angle 
of the present delta wing) or in capturing a shock wave (be- 
cause of the relatively coarse grid distribution in the stream- 
wise direction). The present grid is found to give a reasonable 
agreement with the experiment because the leading-edge vor- 
tex is formed at relatively low angles of attack due to the 
sharp leading edge. 

The computed results are shown here for flow conditions 
at Af, = 0.9 and a = 4 deg. The Reynolds number based on 
the root chord is about 17 million. Figure 3 shows the com- 
parison of computed steady pressures using the modified and 
unmodified Baldwin-Lomax turbulence models with the ex- 
perimental data at 34, 54, and 68 % semispan sections. Since 
the modification accounts for the leading-edge separation, the 
pressure distribution gives a suction peak due to the leading- 
edge vortex. In contrast, the unmodified model smears out 
the vortex so that the suction peak due to the vortex disap- 
pears. Thus, the computed results with the modification of 
the turbulence model successfully capture the flow structures 
shown by the experimental data. 

In the plots, the computed results show good agreement 
with the experimental data at inboard sections. On the other 
hand, the comparisons at the outboard section show that the 
peaks near the leading edge in the computed profiles are 
located farther downstream than in the experimental data. In 
other words, the computed results predict the location of the 
leading-edge vortex more inboard than the experiment. (As 
a result of the discrepancy in the leading-edge region, the 
downstream upper-surface pressure shows a minor disagree- 
ment, too.) 

Since a fine-grid calculation using 151 x 41 x 41 points 
did not improve the comparison at the leading edge, the source 
of the discrepancy may be in the computational model rather 
than the grid resolution. One of the possible sources is the 
wall effect of the wind tunnel. The present computation as- 
sumes a freestream at the far field. Another source is use of 
a boundary-layer transition strip in the experiment. Even though 
the modified Baldwin-Lomax model is used here, the com- 
putation assumes a fully turbulent flow. The other possible 
source is the geometry definition at the tip. The report 9 does 
not have enough information to represent a detailed wing tip 
shape. In addition, the grid represents a clean wing, whereas 
the experimental model has some surface irregularities be- 
cause of control surfaces both at the leading edge and at the 
trailing edge. 
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In Ref. 19, four steady-state cases are computed for flow 
conditions at A/ x = 0.88 and 0.9 with angles of attack a = 

3 and 4 deg. The experimental data do not show the presence 
of a shock wave at A/ x = 0.88. With an increase of the Mach 
number from 0.88 to 0.9, a shock wave is formed on the upper 
surface of the wing. On the other hand, an increase in the 
angle of attack from 3 to 4 deg at a fixed Mach number 
primarily affects the strength of the leading-edge vortex. The 
computed results successfully represent the effects of the dif- 
ferent Mach numbers and the different angles of attack. 

Rigid Pitching Motion 

The unsteady data are given for the case when the rigid 
wing is oscillating in a pitching mode, a(t) = a m - ds\n((ot) y 
about an axis at 65.22% root chord, where cj is the circular 
frequency in radians per second. The test cases consider four 
flow conditions at = 0.88 and 0.9 with mean angles of 
attack a m = 3 and 4 deg, a pitch amplitude a = 0.5 deg, 
and a frequency of 8 Hz that corresponds to a reduced fre- 
quency of k — 0.6 ( k = a jc/I/oc where c is the root chord). 

Unsteady computations are started from the corresponding 
steady-state solutions. The number of time steps per cycle of 
3600 was chosen from the numerical experiments to assure 
the time accuracy (the typical time-step size was about 3.3 x 
10 ~ 3 ). The convergence of the unsteady computations to a 
periodic flow is verified by comparing the results between 
cycles. The third-cycle results are shown in the following. The 
numerical transient is confirmed to disappear within two cycles. 

Figure 4 shows the comparison of computed unsteady pres- 
sures using the modified and unmodified Baldwin-Lomax tur- 
bulence models with the experimental data at M , = 0.9 and 
a m = 4 deg, corresponding to Fig. 3. The plots show the 
comparison of the magnitude and phase angle between the 
computed and measured unsteady upper surface pressure 
coefficients of the wing at 34, 54, and 68% semispan sections. 
Both the leading-edge vortex and the shock wave produce a 
peak in magnitude and a jump in phase angle. Since the 
unmodified model gave almost two orders of magnitude higher 
turbulent viscosity at the leading edge, the solution became 
highly dissipative and thus did not show any large changes in 
unsteady pressures. (In Fig. 4a, the modified turbulence model 
predicts a larger phase change than the experiment near xlc 
= 0.15. This is partly due to the conversion of unsteady 
pressures from real and imaginary to magnitude and phase 
angle.) The improvements due to the modification of the 


turbulence model are seen in both magnitude and phase angle 
where the leading-edge vortex exists. Consistent to the steady 
pressures, the peaks near the leading edge in the computed 
profiles are located more downstream than the experimental 
data at the outboard sections. Overall, the numerical results 
show fairly good agreement with the experimental data. The 
present grid leads to reasonable resolution for the present 
unsteady flowfields, even though the grid is fairly coarse. 

Throughout the four test cases presented in Ref. 19, the 
modified turbulence model predicted higher peaks in mag- 
nitude and larger changes in the phase angle at the leading- 
edge vortex and thus agreed with the experiment better than 
the unmodified model. The effect of Mach numbers on the 
shock wave and the effect of angles of attack on the vortex 
were consistent to the steady-state results. 

Rigid and Flexible Ramp Motions 

In maneuvering, aircraft often undergo rapid ramp mo- 
tions. During such motions, flow unsteadiness and wing flex- 
ibility play important roles. In this section, the applicability 
of the present development to computing such flowfields is 
demonstrated. 

Computations are performed for rigid and flexible wings in 
ramp motion. Structural properties of the wing were selected 
to represent a typical fighter wing. Figure 5 shows the mode 
shapes and the frequencies of the first four normal modes for 
the clipped delta wing used in the following computations. 
The dynamic pressure is set to be 1.0 psi. Test cases consider 
4- and 10-deg ramp motions from 0-deg angle of attack for 
both rigid and flexible wings. 

4-Deg Ramp Motion 

Figure 6 shows the comparisons of the sectional lift re- 
sponses between the rigid and flexible wings at = 0.9 and 
Re c = 15 x 10 6 for the wing ramping up to 4 deg with a pitch 
rate of A = 0.04. The pitch rate A is defined as ac/U *. The 
variation of the effective angle of attack including both the 
ramp angle and the flexible angle of attack is also shown for 
the flexible-wing case. The data are plotted at 34, 54, 68, and 
90% semispan sections. The unsteady computations are started 
from the converged steady-state solution at 0-deg angle of 
attack. 

In the rigid-wing case, the lift responses at the inboard 
sections settle down quickly after the ramp motion stops, and 
the flow approaches the steady-state values. (Thus the com- 
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Fig. 12 Instantaneous density contour plots at 68% semispan section in 
10 6 , A = 0.04. 

putation was stopped at t - 0.5 s.) At the 90% section, how- 
ever, the lift continues to increase for a short period after the 
ramp motion stops. This corresponds to the movement of the 
leading-edge vortex as indicated in the corresponding pressure 
history shown in Fig. 7. In Fig. 7, the pressure distributions 
are plotted every 100 time steps. The ramp motion ends at 
600 time steps ( t = 0.04 s), which corresponds to the sixth 
plot from the bottom. At the 90% section in Fig. 7d, a leading- 
edge vortex is formed and lifts off from the wing surface. This 
results in the increase and then a decrease of the sectional lift 
shown in Fig. 6. At the other inboard sections in Fig. 7, the 
pressure distributions become similar to those in Fig. 3, i.e., 
the leading-edge vortex remains without interacting with the 
shock wave. Thus, the corresponding sectional lift stays nearly 
steady with minor fluctuations due to the movement of the 
vortex, as seen in Fig. 7a. 

In contrast to the rigid-wing case. Fig. 6 shows that the 
sectional lift responses of the flexible wing are oscillatory. 




10-deg ramp motion (0.08 < t < 0.16 s), M x = 0.90, Re c = 15.0 x 

The primary oscillation is about 6 Hz, which corresponds to 
the frequency of the first mode. As the wing moves down- 
ward, the lift increases, and vice versa. The overall oscillation 
is damping. The first mode response is damping, but the sec- 
ond mode response stays oscillatory. 19 The primary effect of 
the flexibility is the reduction of the lift at most of the span- 
wise sections due to the reduced angles of attack. Figure 8 
shows the corresponding pressure history plots. The plots in 
Figs. 8c and 8d show the shedding of the leading-edge vortex 
in comparison to the plots of the rigid wing in Figs. 7c and 
7d. In addition, Fig. 6 shows that the 90% sectional lift re- 
sponse of the flexible wing has the low-frequency primary 
oscillation perturbed by a high-frequency vortex shedding. 
The low-frequency primary oscillation (6 Hz) corresponds to 
the structural oscillation. The high-frequency perturbation 
corresponds to the shedding of the vortex that can be seen in 
the pressure plots of Fig. 8d. This frequency (about 15 Hz) 
does not excite any structural mode as shown in Fig. 5. 
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Clipped Delta Wing (Flexible) 



Fig. 13 Streamline pattern over the upper surface of flexible wing 
in 10-deg ramp motion, = 0.90, Re c = 15.0 x 10 6 , A = 0.04. a) 
t = 0.10 s, and b) / = 0.16 s. 

Computations (not shown) were also carried out for the 
flexible-wing ramping from 0 to 3 deg and from 0 to 5 deg at 
the same pitch rate as the 4-deg case. The 4-deg case shows 
the largest high-frequency perturbation. In the 3-deg case, 
the vortex is not strong enough to disturb the lift response. 
The vortex shedding is found only at the 90% section. In the 
5-deg case, the vortex lifts off from the wing surface so that 
the structural oscillation does not cause the perturbation seen 
at the lower angles of attack. The 10-deg case discussed in 
the next section does not show the perturbation either. 

A reduction of the local angle of attack due to the flexibility 
of the wing results in a delay of the lift increase for a short 
period after the ramp motion stops (about 0.04 < t < 0.11 in 
Fig. 6). The wing is deformed upward in bending and leading- 
edge down in twisting. When the ramp motion stops, the wing 
is still deforming, which also gives a dynamic effect. Thus the 
local angle of attack relative to the wing section decreases 
toward the wing tip and the leading-edge vortex appears weaker 
in the flexible-wing case. This leads to the delay of the lift 
increase in the flexible-wing case. 

10-Deg Ramp Motion 

The sectional lift responses for the 10-deg ramp motion at 
several spanwise sections are shown in Fig. 9. The compu- 
tations are again started from the steady-state solution at 0- 
deg angle of attack. The unsteady increase of the lift is ob- 
served more widely in both rigid and flexible cases than the 
4-deg case. The sectional lift at the 90% section indicates a 
stall before reaching 10-deg angle of attack for both rigid and 
flexible wings. Instead, the plot does not have any significant 
perturbations. The flexible wing gives lower lift because of 
the deformation of the wing similar to the 4-deg case. After 
the initial lift increase, the lift oscillates due to the structural 
oscillation. Again, the first mode response is damping, but 
the second mode response stays oscillatory. 19 Figure 10 shows 


the magnified deformation of the wing at 1600 time steps (the 
ramp motion ends at 1500 time steps: t = 0.10 s). The actual 
displacement of the leading edge at the wing tip is 1.7% of 
the root chord length. 

Figure 11 shows the corresponding pressure history plots 
of every 100 time steps for the flexible-wing case. In contrast 
to the 4-deg case, the deformation of the wing does not affect 
the flowfield as strongly because the leading-edge vortex lifts 
off from the wing surface at the outboard sections. The flex- 
ibility does not play an important role at the inboard sections, 
because the wing root is fixed. Thus there is no significant 
difference in the responses between the rigid and flexible 
wings. At the most inboard section in Fig. 11, the pressure 
distributions show no interaction of the leading-edge vortex 
with the shock wave. At the 54% section, both vortex and 
shock wave develop, then both disappear. At the 68% section, 
a similar but more rapid change occurs. This rapid reduction 
of the lift indicates a vortex breakdown. 

To see the interaction of the leading-edge vortex with the 
shock wave, the density contours at the 68% section are plot- 
ted every 200 time steps from 1200 to 2400 time steps (0.08 
< t < 0.16) in Fig. 12. First, there is no interaction between 
the vortex and the shock wave. As the vortex develops, it 
moves toward the trailing edge, that is, toward the shock 
wave. When the vortex starts to interact with the shock wave, 
the shock wave starts to ride on the shear layer and to form 
a lambda-type shock wave. At this point, the shock wave 
disappears from the surface pressure plots. As the front shock 
grows, the flow separation grows and the vortex core bursts 
quickly. Simultaneously, the rear shock weakens. Finally, the 
fully separated flow is observed. The corresponding contour 
plots of the negative u component in Ref. 19 show that a 
negative u region appears as the vortex is deformed by the 
strong front shock at 2000 time steps ( t — 0.13 s). As the 
vortex core diffuses, the reverse flow region grows. The shock 
wave plays an essential role in the process of this breakdown. 

Figure 13 shows plots of the streamline pattern at 1500 and 
2400 time steps ( t = 0.10, 0.16 s, respectively). At 1500 time 
steps, when the ramp motion has just ended, a leading-edge 
separation is formed over the entire span. Although it is not 
clear from the plots, the vortex starts bursting near the tip as 
indicated in Fig. 1 Id. In contrast, at 2400 time steps, a bubble- 
type breakdown is clearly observed in the middle of the span. 
The breakdown grows toward the upstream slowly. Then the 
flow reaches the nearly steady state. 

Although no corresponding experiment was performed for 
the present wing at this angle of attack, Ref. 20 reported a 
vortex breakdown about a similar wing at similar flow con- 
ditions, including a comparison with the experiment. 

Concluding Remarks 

In this paper, a computational procedure for computing the 
unsteady transonic flows associated with the leading-edge vor- 
tex on a clipped delta wing, including flexibility, has been 
presented. The procedure is based on a time-accurate com- 
putational method combined with the use of aeroelastically 
adaptive dynamic grids. The flow is modeled using the Navier- 
Stokes equations. The flow equations are coupled with the 
structural equations to account for the flexibility. The nu- 
merical procedure has been verified through the comparisons 
with the experiment for the unsteady pitching cases on the 
rigid clipped delta wing. The main flow structures are suc- 
cessfully captured. 

The ramp motion cases have demonstrated the effects of 
unsteadiness of the flow field and flexibility of the wing. The 
primary effect of the flexibility is the reduction of the lift due 
to the deformation of the wing. Interaction of the leading- 
edge vortex with the shock wave has significant effects on the 
wing responses. For the 4-deg ramp motion, the vortex shed- 
ding occurs at the wing tip due to the flexibility. For the 10- 
deg ramp motion, a possible vortex breakdown is observed. 
The inviscid interaction with the shock wave plays an essential 
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role in the process of the breakdown observed in the present 
calculation. 
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Abstract 

Unsteady Navier-Stokes computations have been 
performed for simulating transonic flows over wings with 
oscillating control surfaces using a locally moving grid 
and a stationary-mismatched zoning scheme. An F-5 wing 
and a clipped delta wing are chosen for the present study. 
The computed unsteady pressures and the response 
characteristics to the control surface motions are compared 
with experimental data. The results successfully predict 
main features of the unsteady pressure profiles, such as the 
double peaks at the shock wave and at the hinge line. 

Introduction 

Aircraft are often subject to aerodynamic oscillation, 
especially in the transonic regime, because of flow 
nonlinearities and the presence of the moving shock 
waves. In this unsteady aerodynamics environment, 
aircraft rely heavily on active controls for safe and steady 
flight operation. Active control is also needed for the 
suppression of structural flutter and the reduction of 
structural weight to achieve stable flight conditions. 

The influence of control surfaces on both 
aerodynamics and aeroelastic performance of a wing is 
more pronounced in the transonic regime. These 
influences can be constructively used to improve the wing 
performance through proper maneuvering of the active 
control surfaces. Active control technology relies on 
accurate predictions of unsteady aerodynamics and 
aeroelastic performance of a wing. Since the experimental 
evaluation of the effect of a control surface on the wing 
performance would involve considerable cost and the risk 
of structural damage in a wind tunnel, it is necessary to 
initiate the investigation through theoretical analyses. 

The theoretical analysis of transonic flows is 
complicated by the presence of mixed subsonic and 
supersonic regions within the flow field. For an unsteady 
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flow field, such as that surrounding a control surface, 
additional considerations are needed to treat moving shock 
waves of varying strength and subsequent flow separation 
induced by the shock-wave/boundary-layer interactions. 
For the case of the control surface in which viscous effects 
dominate, computation based on the unsteady Navier- 
Stokes equations is needed. 

The physics of unsteady transonic flow around a 
control surface has been simulated at various levels of 
inviscid and viscous approximations using small 
disturbance theory 1 * 2 and, for limited two-dimensional 
cases, the unsteady Navier-Stokes equations. 3 * 4 The 
purpose of this study is to explore the capability of three- 
dimensional Navier-Stokes simulation for the unsteady 
flow field surrounding a wing with an oscillating control 
surface. 

The present investigation is initiated in conjunction 
with a recently developed code, ENSAERO, which is 
capable of computing aeroelastic responses by 
simultaneously integrating the Euler/Navier-Stokes 
equations and the modal structural equations of motion 
using aeroelastically adaptive dynamic grids. 5 ' 8 The code 
has been applied to transonic flows from small to 
moderately large angles of attack for fighter wings 
undergoing unsteady motions. In this paper, the geometric 
capability of the code is extended to simulate unsteady 
flows over a rigid wing with an oscillating trailing-edge 
flap. 

To model an oscillating control surface efficiently, an 
algebraic grid generation technique is incorporated into the 
code. The grid moves every time step to follow the 
deflection of the control surface. Small deflections are 
handled using a sheared single grid, and large deflections 
are handled using a zonal technique. 9 

In this paper, the first test case considers transonic 
flows over an F-5 wing with an oscillating inboard 
control surface. The same case was simulated using small 
disturbance theory in Ref. 2. The second case considers 
transonic vortical flows over a clipped delta wing. 
Unsteady Navier-Stokes computations for the clean wing 
were reported in Ref. 8. The mismatched zonal cases are 
demonstrated for this case. 
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Numerical Method 


Governing Equations and Discretization 

The nondimensionalized Reynolds-averaged thin-layer 
Navier-Stokes equations used in this study can be written 
in conservation-law form in a generalized body- 
conforming curvilinear coordinate system for three 
dimensions as follows: 


See Ref. 13 for more details. Because control surfaces 
oscillating at a small amplitude were considered in this 
paper, the FD and FV formulas did not produce any 
significant difference in the results. 


3 t Q + 3 eE + 3„F + drG - 9^G V (1 ) 

Re 

where t - t, 1 = |(x,y,z,t), n = T|(x,y,z,t), and C - 
(^(x.y.z.t). In the present paper, the | and n directions are 
along the streamwise and spanwise directions, 
respectively, of a wing. The viscous derivatives associated 
with these directions are dropped. The C direction is 
normal to the wing surface, and thus the viscous 

derivatives are retained. 

/\ ^ 
In Eq. (1), Q is the vector of conserved quantities, E, 

/N ✓N — V 

F and G are the inviscid flux vectors, and G is the thin- 
layer viscous flux vector. The viscosity coefficient in 

G is computed as the sum of the laminar and turbulent 
viscosity coefficients where the laminar viscosity is taken 
from the freestream laminar viscosity, assumed to be 
constant for transonic flows. As an option, Sutherland's 
law can be used to calculate the laminar viscosity. The 
turbulent viscosity is evaluated by the Baldwin-Lomax 
algebraic eddy-viscosity model. 10 Since the flow field to 
be considered in this paper contains leading-edge 
separation, it is important to apply a modification to the 
turbulence model originally developed for crossflow-type 
separation. 11 

In this study, a conventional finite-difference (FD) 
grid system is used. However, for the moving grid case, it 
is important to treat the time metrics correctly. The FD 
formulation does not give freestream-capturing time 
metrics. 1213 The FD time metrics are related to the grid 
velocity, for example. 

It - “ lx*t - lyYt - IzZt (2) 

For rigid motion of the grid, the freestream can be 
preserved by applying a freestream subtraction technique 
to the time metric terms. 

In the FV formulation, the right-hand side of Eq. (2) 
is expressed as a product of a surface vector and a grid 
velocity. The time integration of the product from time 
level n to n+1 represents a volume swept by the surface 

during one time step. By computing this volume, 
the freestream-capturing time metrics can be obtained as 


Upwind Algorithm 

Several numerical schemes have been developed to 
solve Eq. (1). The present code has two different schemes: 
the central-difference and streamwise upwind schemes. In 
this work, a streamwise upwind algorithm is applied to 
compute the inviscid cell-interface fluxes. The streamwise 
upwind algorithm has recently been developed and applied 
to steady and unsteady problems of transonic flows over 
wings, including flexibility of the wing. 7 - 8 A third-order 
evaluation is used for the inviscid term. A second-order 
central -difference evaluation is applied to the viscous term. 
The complete algorithm can be found in Ref. 7. 

An implicit method is used for the time integration 
because it is more suitable for expensive unsteady viscous 
calculations. The method chosen for the upwind scheme is 
the LU-ADI (lower-upper factored, alternating direction 
implicit) method. 14 The method is first-order accurate in 
time and requires only scalar bidiagonal matrix inversions. 
See Refs. 14 and 15 for additional details. 

Specific code performance information for the current 
study is given as follows. All results were computed on 
either a CRAY-YMP or CRAY-2 computer at NASA 
Ames Research Center. The performance of the upwind 
version of ENSAERO is 170 MFLOPS and 19 psec per 
iteration per grid point on a single CRAY-YMP 
processor. 

Control Surface Grid 

In this paper, a single C-H topology grid for a wing 
with a control surface is used. In each C grid section, a 
two-dimensional algebraic grid generation technique is 
applied to each airfoil section. 16 At both ends of the 
control surface, a small gap is introduced (Fig. 1). This 
region is used to shear the grid when the control surface 
oscillates. Although the gap is introduced to simplify the 
calculations, its effect can be minimized by clustering the 
grid in this region. 

The C grid around a deflected control surface can be 
obtained in two ways. One is to shear every grid line 
normal to the control surface with the local deflection. Ax 
and Az. The other is to regenerate the entire C grid with 
the control surface deflected at every time step. The 
computational overhead of the latter approach is only 7% 
because of the efficient algebraic grid generation scheme. 
However, the computed surface pressures did not show any 
differences between the two methods. Thus, the results 
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presented in this paper were computed using the shearing 
technique. 

Zonal Grid Capability 

The present code inherits the zonal grid capability 
developed for the Transonic Navier-Stokes (TNS) code to 
handle complicated geometries, such as complete aircraft 
configurations. 9 Although a problem with the extension 
of the zonal method to unsteady flows was reported in a 
previous paper, 17 a code error was found and has been 
corrected. In the present study, the corrected code is used. 

In this paper, a mismatched zoning scheme is 
introduced to accommodate iarge mean flap deflections, 
while a shearing-grid technique is used to model an 
oscillating flap about the mean deflection. Zonal interfaces 
are located at both ends of a control surface including the 
gap (Fig. 2). As the control surface deflects, the grids 
become mismatched. However, if the zonal interfaces 
move, it requires expensive computation to find 
interpolation coefficients on the mismatched zones at 
every time step. To maintain the efficiency of the single- 
grid computation, the zonal interfaces should remain 
stationary when the control surface oscillates with a small 
amplitude. Instead, the control surface grid shears at the 
gap region similar to the single-grid case. To transfer the 
flow information from one grid to another, bilinear 
interpolation is used here (see the nonconservative 
interpolation in Ref. 18), because the present zonal 
interfaces are coplanar. Although the present interpolation 
is explicit and nonconservative, the error can be ignored 
when practical time-step sizes for Navier-Stokes 
computations are used. 

A different zonal approach is reported in Ref. 19 to 
treat a stationary, deflected flap. Its extension to 
oscillating control surfaces will be investigated in the near 
future. 

Results 

F-5 Wing 

The first test considers unsteady viscous flows over 
an F-5 wing with an oscillating inboard control surface. 
This wing has an aspect ratio of 2.98, a taper ratio of 0.31 
and a leading edge sweep angle of 31.92 deg. 
Computations were made using coarse and fine grids 
containing 151 x 41 x 34 points and 201 x 61 x 34 
points, respectively. The wing planform is given in Fig. 

3. It should be noted that the present C-H grid does not 
have enough resolution for the faired wing tip of the 
experimental model. The control surface is oscillating 
about an axis located at the 82% root chord, and the hinge 
axis is normal to the wing root. The test cases are at a 
Mach number of M « 0.9, where the experimental steady 
and unsteady data are given in Ref. 20. All F-5 wing cases 
are computed at a Reynolds number based on the root 


chord olRe c » 12 x 10 6 . Reference 2 discusses the small 
disturbance results applied to the same test case. 

The response of surface pressure to the control surface 
motion can be represented m terms of real and imaginary 
parts of the first Fourier component of the unsteady 
pressures. In Fig. 4, the coarse- and fine-grid results are 
compared with experimental data at an angle of attack of a 
m ® with the control surface oscillating at a frequency 
of 20 Hz and an amplitude of p - 0.5 deg. This frequency 
corresponds to a reduced frequency of k = 0.28. Results are 
shown for the upper surface pressures at three spanwise 
locations. It is noted that for * 0.9, the steady-state 
solution is shock-free except near the wing tip. The spikes 
in the unsteady pressure distributions around the 50% 
chord indicate the motion of the shock wave due to the 
control surface oscillation. The spikes are also seen at the 
hinge line. In the real part of the unsteady pressures, at the 
inboard sections, the computations predict higher spikes 
for the motion of the shock wave than observed in the 
experiment. This is because the computation assumes a 
plane of symmetry at the root section, while the 
experiment has a solid wall. In contrast, the computations 
predict lower spikes at the hinge line. This is because the 
present grid has constant chordwise distributions and thus 
does not align to the hinge line. In the imaginary part, 
there is a greater discrepancy between the computation and 
the experiment, which is possibly due to the resolution of 
the experimental data as shown later. The coarse- and fine- 
grid results show reasonably good agreement in Fig. 4. 

The unsteady pressure profiles using time-step sizes 
yielding 1 800 and 2400 steps/cycle were compared with 
each other to check the time-step dependency of the coarse* 
grid results. The unsteady results converged at 1800 
steps/cycle and thus this time step was used for the F-5 
wing results shown in this paper. The finite-difference 
time metrics, Eq. (2), and the freestream-capturing time 
metrics, Eq. (3), were checked on the coarse grid as well. 
No difference was found for this small amplitude of 
control surface oscillation. 

Figure 5 shows the unsteady results on the coarse grid 
at 1.5 deg angle of attack with the flap oscillating at a 
frequency of 20 Hz and an amplitude of 0.5 deg at three 
spanwise locations. The real part of the unsteady pressures 
follows the previous observation for Fig. 4. The 
imaginary part shows better agreement at the outboard 
section than Fig. 4. Since the computed results show 
more consistent trends, the disagreement in the imaginary 
part is likely due to the resolution of the experimental 
data. 

Clipped Delta Wing 

The next test case considers a clipped delta wing with 
an oscillating trailing-edge control surface. 21 The wing 
planform is shown in Fig. 6. The wing has a leading-edge 
sweep angle of 50.4 deg and a 6%-thick circular-arc airfoil 
section. At » 0.9 and tx ■ 3 deg, both a leading-edge 
vortex and a shock wave are present on the upper surface 
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of the wing. The present C-H grid contains 15 1 x 44 x 34 
points (see Fig. 1). Since the experiment was conducted 
using a Freon test medium, the ratio of specific heats, y , 
is set to 1.135 in the following computations. In addition, 
the modified Baldwin-Lomax model is used to account for 
the leading-edge vortex. Steady-state and rigid pitching 
calculations of this wing were reported in Ref. 8. 
Although three cases were computed for this oscillating 
control surface with an amplitude of 2.20, 4.38 and 6.65 
deg, the 6.65-deg case is discussed in the following 
because the other two cases show similar trends. 

Figure 7 shows the unsteady pressures with the 
control surface oscillating at a frequency of 8 Hz and an 
amplitude of 6.65 deg at - 0.9, a - 3 deg and Re e - 
17xl0 6 based on the root chord. Results are shown as 
magnitude and phase angle of the upper surface pressure 
responses at three spanwise sections. The magnitude part 
of the unsteady pressures shows significant influence of 
the control surface oscillation. At the 68% spanwise 
section, the plot shows double peaks, similar to the F-5 
wing case. These correspond to the shock wave location 
and the control-surface hinge line. Although the computed 
result shows a much smaller peak due to the shock wave 
at the 68% section, it is consistent with the steady-state 
result where the shock wave appears weaker. 8 This 
discrepancy originates in both the experimental and 
computational modeling. For example, there was a 
misalignment of the control surface on the order of 2 deg 
in the experiment. 21 In contrast, the computed peak at the 
hinge line is much higher than the experiment. This is 
because the grid aligns to the hinge line in this case and 
the experimental data does not have enough resolution 
near the hinge line. The phase angle plots show a 
discrepancy near the leading edge. It should be noted that 
the computation assumes a fully turbulent flow while the 
experiment has a fixed transition at the 8 % chord. 
Overall, the computed results show reasonably good 
agreement with the experiment. 

The test case presented in Fig. 7 was more sensitive 
to time-step sizes. Most of the unsteady pressures on the 
wing surface converged at 3600 steps/cycle (the same 
number as the rigid pitching case in Ref. 8). However, the 
oscillation of the control surface caused a minor 
fluctuation to the leading-edge vortex, and thus, the 
unsteady pressures under the leading-edge vortex required 
5000 steps/cycle to converge. 

Zonal Grid 

The zonal test case is taken from the previous clipped 
delta wing. The control surface is deflected 6 deg and 
undergoes an oscillatory motion with an amplitude of 1 
deg and a frequency of 8 Hz at - 0.9, a - 0 deg and 
Re c - 17x1 0 6 . Since no experimental data is available for 
this particular case, zonal-grid results are compared with 
single-grid results. 

First, the single grid of 151 x 44 x 34 points is 
divided into three zones at the gaps at both ends of the 


control surface. Each zone contains 151 x 13 x 34 , 151 x 
15 x 34, and 151 x 20 x 34 points from inboard to 
outboard (see Fig. 2). Special attention is given to the 
generation of smooth zonal interfaces because the gap 
region often causes skewed grid distributions. 

Figure 8 shows the comparison of steady-state results 
between the single- and zonal-grid solutions. The control 
surface is deflected down 6 deg and held stationary. 
Asymmetric pressure distributions indicate the effects of 
the deflected control surface. These effects can even be 
seen at the inboard sections. A comparison of unsteady 
pressures between the single- and zonal-grid results is also 
given in Fig.^9. The control surface oscillates with an 
amplitude of p - 1 deg about a mean deflection of p m - 6 
deg. The unsteady pressure profiles show double peaks at 
the shock wave and the hinge line, although the shock 
wave is not seen in Fig. 8. In both figures, single- and 
zonal-grid computations agree' well, in spite of the 
mismatched zoning. It takes about 9 sec on a CRAY- 
YMP single processor to set up the interpolation 
coefficients for all of the zonal interfaces. This 
corresponds to approximately 2 time steps of the Navier- 
Stokes calculation. After the coefficients are computed, 
the unsteady computations simply use the stored 
coefficients because the zonal interfaces remain stationary. 
This is one way to take an advantage of the flexibility of 
the zonal-grid technique along with the efficiency of the 
single-grid computation. 

Conclusions 

Unsteady Navier-Stokes computations have been 
performed for transonic flows over wings with oscillating 
control surfaces. To use the existing framework of the 
finite-difference method, the time metric terms are 
calculated by using the finite-volume concept. The use of 
locally moving grids has made it possible to efficiently 
simulate an oscillating control surface. The introduction 
of the stationary-mismatched zoning scheme has enhanced 
the flexibility of the code. 

Comparisons of unsteady pressures with experimental 
data show reasonably good agreement, although the grids 
used in this study are relatively coarse. The computed 
results successfully predict main features of the unsteady 
pressure profiles, such as the double peaks at the shock 
wave and at the hinge line. In the future, an unsteady 
mismatched zoning scheme will be investigated. A control 
law will also be implemented into the code to compute 
aeroelastic performance of a wing. 
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Fig. 1 A single grid for a wing with a control surface. 



Wing Trailing Edge 


Control Trailing Edge 


Gap Region 
(Enlarged Scale) 


Gap Region 
(Enlarged Scale) 


Wing Trailing Edge 


Crossflow View at the Trailing Edge 


Fig. 2 Zonal grids for a wing with a control surface. 
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Fig. 3 Planform of F-5 wing with inboard trailing-edge flap. 


F-5 wing Fine grid 

Moo = 0.9, Re = 12 :< 10^ (201 x 61 x 34) 

a = 0°, 5 = 0.5° Coarse grid 

k = 0.278 0 51 x 41 x 34 ) 

□ Experiment, Persoon et ai. 


F-5 wing Computation 

Moo = 0.9, Re = 12 x 10 6 D Experiment, Persoon et al. 

a =1.5°, 5 = 0.5° 
k = 0.28 

Grid: 151 x 41 x 34 


Semispan 

84% 


Semispan 

84% 




Fig. 4 Comparison of computed unsteady pressures Fig. 5 Comparison of computed unsteady pressures with 
between the coarse- and fine-grid solutions on the upper experimental data on the upper surface of the F-5 wing, 
surface of the F-5 wing. 
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Fig. 6 Planform of clipped delta wing with trailing-edge 
flap. 


Clipped delta wing Computation 

Moo = 0.9, Re = 17 x 10 ® □ Experiment, Hess et al. 

a = 3°, 5 = 6.65° 
k = 0.588 

Grid: 151 x 41 x 34 

Semispan 

68% 



Fig. 7 Comparison of computed unsteady pressures 
between the use of freestream-capturing and finite- 
difference time metrics on the upper surface of the clipped 
delta wing. 


311 


Steady C, 


Clipped delta wing Single-grid computation 

M m = 0.9, Re = 17 x 10 6 (151 x 44 x 34) 

a = 0° 8.- - 6° Zonal-grid computation 

’ ,lx “ (151x13x34, 

151 x 15 x 34, 

151 x 20 x 34) 



0 -2 .4 .6 .8 i.o 

X/C 


Fig. 8 Comparison of computed steady pressures between 
the single and zonal grids over the clipped delta wing. 
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Abstract 

A new unsteady zoning method called ‘virtual zones’ has been developed for application 
to an unsteady finite difference Navier-Stokes code. An existing interpolation method has 
been extensively modified to bring the run times for the interpolation procedure down to the 
same level as for the flow solver. Unsteady Navier-Stokes computations have been performed 
for transonic flow over a clipped delta wing with an oscillating control surface. The computed 
unsteady pressure and response characteristics of the control surface motion compare well 
with experimental data. 


Introduction 

Present civilian transport aircraft as well as highly maneuverable fighter aircraft are 
often subject to unsteady aerodynamics. In this unsteady environment aircraft designers 
utilize active controls to achieve controllability and safety of the aircraft. Active control can 
also be used to suppress transonic flutter characteristic of high aspect ratio transport wings 
and thus reduce the structural weight to achieve more efficient flight conditions. 

In the transonic flow regime active controls have a pronounced effect on the aerodynamic 
and aeroelastic performance of a wing. This effect can be used to improve the airplane 
performance by the proper design of the active control surfaces. To do this successfully 
requires the accurate prediction of the aerodynamic and aeroelastic performance of a wing. 
Experimental prediction of unsteady aerodynamic and aeroelastic performance is costly and 
time consuming; numerical simulation of the unsteady Navier-Stokes equations is a much 
more cost effective alternative for predicting the performance of an active control surface. 

The physics of unsteady transonic flow around a control surface has been simulated with 
small disturbance theory [1,2]. Unsteady Navier-Stokes simulations have been done in two 
dimensions [3,4]. A recent study [5] explored a three dimensional simulation of the unsteady 
thin-layer Navier-Stokes of the flow field surrounding a wing with an forced oscillating control 
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surface. In that study an unsteady Navier-Stokes code, ENSAERO, was extended to simulate 
unsteady flows over a rigid wing with an oscillating trailing-edge flap. 

In the previous study an algebraic grid generation technique was incorporated into the 
code. The grid moved at every time step to follow the deflection of the flap. The small 
unsteady deflections were handled using a sheared single mesh. The large stationary deflection 
were handled using a zonal method [6]. The use of the single sheared grid did not permit the 
exact simulation of the unsteady flap-wing geometry. A gap had to be introduced between 
the ends of the flap and wing to allow sufficient space for the moving sheared mesh. The 
gap compromised the numerical simulation of the oscillating control surface flow field. The 
purpose of the present study is to rectify that compromise through the use of a new zoning 
technique called “virtual” zones. 

The numerical simulation of the unsteady Navier-Stokes equations about complex and 
realistic aerodynamic configurations requires the use of zonal methods. In this method the 
overall flow field domain is subdivided into smaller blocks or zones. In each of these zones 
the flow field is solved independently of the other zones. The boundary data for each zone 
is provided by the neighboring zones. A major difficulty of the zonal methods applied to 
oscillating control surfaces has been how to account for the variable exposure of the ends of 
the control surfaces to the flow field. 

The virtual zoning method, first implemented in a multizone finite volume code, CNSFV, 
[7], has been modified for application to the unsteady finite difference code, ENSAERO. For 
a finite difference application, virtual zones are essentially two-dimensional zones one cell 
thick. The main purpose of these zones is to convert, for example, a solid wall boundary 
condition into an interface condition. The interface conditions are required for the interzonal 
communication. In a multi-zonal code the virtual zones act like real zones as far as boundary 
and interface conditions are concerned, however, no flow field computations are done within 
these zones. Hence the name ‘virtual’ zone is appropriate. 

In addition to the introduction of the virtual zones, it is necessary to speed up the process 
of determining the interpolation coefficients required for the interzonal communication if the 
unsteady Navier-Stokes simulation is to be practical. 

The present study considers the transonic vortical flow over a clipped delta wing. A view 
of the wing and the control surface is shown in Fig. 1. Unsteady Navier-Stokes computations 
for the clean wing were reported in Ref. 8. The forced oscillating control surface computation 
with the single zone sheared mesh were presented in Ref. 5. 

Numerical Method 


Governing Equations and Discretization 

The governing equations are the Reynolds- averaged thin-layer Navier-Stokes equations. 
The laminar viscosity is taken from the freestream laminar viscosity and is assumed to be 
constant for the transonic flow considered in this study. The turbulent viscosity is obtained 
with the Baldwin- Lomax algebraic eddy viscosity model [9] with the Degani-Schiff modifica- 
tion [10] to properly handle the leading edge separation as well as the control flap vortical 
flow. 

The numerical algorithm, the time dependent metrics of the curvilinear coordinate sys- 
tem, and the performance characteristics of the ENSAERO code have been described previ- 
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ously and will not be repeated here. The interested reader is referred to Ref. 5. 

Control Surface Grid and Zones 

The primary focus of the present study is to demonstrate the feasibility of using dynamic 
zones for the oscillating control surface case and not so much as minimizing the cpu run times 
for the case to be presented. Hence the flow domain was split into only three real zones. Each 
of the zones consists of a C-H topology. The two zonal boundaries were placed at the span 
stations located at the ends of the control surface. Four additional zones (called virtual zones 
and described below) were placed in the two cuts separating the flap and the wing. One pair 
of the virtual zones remains fixed with the wing and the other pair is fixed with the flap 
and moves with the flap during the control surface motion. Figure 2 shows the seven zones 
involved. 

The C grid around a deflected control surface can be obtained in two ways. One is to 
shear every grid line normal to the control surface with the local deflection. The other is 
to regenerate the entire C grid with the control surface deflected at every time step with an 
algebraic method. A previous study [5], showed that the computed surface pressures did not 
show any differences between the two methods. Therefore, for this study, the grids around 
the control flap were regenerated with the simpler shearing method. 

Virtual Zone 

The zoning capability of the CNSFV code [6], allowed the possibility of a single face of 
a zone to interact with several other zones. The procedure of determining the interpolation 
coefficients is automatic in that no additional information is required other than identify- 
ing the faces that are in contact with each other. To further extend the flexibility of the 
zoning method for the case of control surface aerodynamics, the idea of "virtual” zones was 
introduced in Ref. 7. 

Virtual zones are zones of zero thickness (for a finite volume formulation) which serve to 
transfer solid wall (or other) boundary conditions to an interface condition. Thus multiple 
boundary conditions can be imposed on a block face with the same flexibility as an interface 
condition. Virtual zones also decouple the process of volume grid zoning from the surface grid 
patches' which define the aerodynamic configuration under study. Surface grid patches are 
required in order to impose the proper boundary conditions. The volume grid zones should 
be set up to obtain the proper mesh qualities required for numerical accuracy. Another 
advantage of the decoupling is that much fewer zones are now needed, thus easing the effort 
and time required to generate the grids about complex and realistic configurations. 

Finally the virtual zones also allow the zonal boundaries to cut through the configuration 
surfaces, which is an important property for control surfaces. The region of the configuration 
that intersects the zonal face is covered with a virtual zone to convert that region into another 
interface condition. Once a zone has been defined along with its associated virtual zones, its 
definition is complete and is not influenced by any of its neighboring zones. In other words 
the zone communicates with the outside world only through the interface conditions. Thus a 
particular zone can be altered or substituted with another zone without any need to redefine 
the interface conditions of the other zones. For example, a zonal grid can be set up for a 
wing with control flaps with one zone for each (say, undeflected) flap. For the deflected flap 
case, only the flap zone needs to be replaced with a zone containing a deflected flap. The 
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boundary and interface conditons of all other zones remain unchanged, even though they 
may now have a solid surface exposed to the flow field, e.g. the edge of the exposed end of 
the wing and flap. 

The orginal zoning capability of the ENSAERO code was extended by including the 
above capability of multiple interface conditions on a single block face. Since the code is a 
finite difference code the zones required an overlap at the boudaries of the zones to allow for 
the proper interblock communication (ie. interfacing). In this study a one cell overlap was 
chosen. For this kind of zoning the virtual zones of zero thickness used for the finite volume 
formulation is not appropriate. Instead the thickness of the virtual zones had to be expanded 
to include the extend of the overlap of the zones. In other words the virtual zones for the 
present formulation are now one cell thick. There is a slight mismatch of a half cell thickness 
between the location of the actual solid wall and the location where the virtual zones applies 
the solid wall boundary condition. This mismatch does not occur with the finite volume 
formulation. However in the present case the slight mismatch had no discernable influence 
on the overall flow field, especially in the case where the ends of the flaps and wings are 
treated viscously, rather than with the inviscid tangency condition. 

An example of the virtual zones required for the present case is shown in Figs. 3-5 for 
the inboard end of the control surface. The two virtual zones slide through each other with 
the control surface motion. As can be seen in the figures different amounts of the virtual 
zones are then exposed to, or in contact with, the real zones surrounding the wing and the 
flap. The virtual zones transfer the solid wall boundary condition to an interface condition 
and thus allows for the automatic inclusion of the variable exposure of the ends of the flap 
and wing to the flow field. The area where the wing and flap virtual zones overlap represents 
the unexposed portions of the flap and wing. Since the flow field is not updated in the 
virtual zones by the flow solver, nothing happens in the virtual zone overlap region nor does 
it influence the rest of the flow field. 

Zonal Interface Interpolation 

The original interpolation procedure used in CNSFV and ENSAERO was based on global 
area search. Even though it was vectorized, it was still much too slow for a dynamic interpo- 
lation procedure where the interpolation coefficients have to be updated at every time step. 
By replacing the area search with a procedure based on a polygon clipping algorithm, the 
search time to find the interpolants was reduced by two orders of magnitude. This improve- 
ment brought the cpu run times for determining the interpolation coefficients down to the 
same level as required for the flow solver (about 4 cpu seconds per iteration on the Cray 
YMP for the case discussed in the next section). More complete details will be presented in 
the full paper. There were some other techniques used that also helped speed up the process 
and these will also be discussed in the full paper. 

Results 

The test case considered in the present study is a clipped delta wing with an oscillating 
trailing edge control surface [11]. The wing panform is shown in Fig. 1. The wing has a 
leading edge sweep angle of 50.4 deg and a 6% thick circular arc airfoil section. At = 0.9 
and a = 3 deg, both a leading edge vortex and a shock wave are present on the upper surface 
of the wing. The C-H grids of the three zones consist of 151 x 13 x 34, 151 x 15 x 34, and 
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151 x 20 x 34 points from the inboard to the outboard zones as shown in Fig. 2. Since the 
experiment was conducted using a Freon test medium, the ratio of specific heats, 7, is set 
to 1.135 in the present computations. As stated before, the modified Baldwin- Lomax model 
is used to account for the leading edge and control surface vortices. Steady state and rigid 
pitching calculations of this wing were reported in Ref. 8. 

Figure 6 shows the unsteady pressures with the control surface oscillating at a frequency 
of 8 Hz and an amplitude of 6.65 deg at M = 0.9, a = 3 deg and Re c = 17 x 10 6 based 
on the root chord. The results are shown as the amplitude and phase angle of the upper 
surface at the three span stations as indicated on the figure. In general the agreement with 
the experimental results is good. Because the accuracy of experiment had its own limitations, 
the virtual zone results are also compared with the single grid results of Ref. 5. As shown 
there is quite a discrepancy between the virtual zone results and the single grid results (the 
151 x 44 x 34 grid), especially in the amplitudes at the center of the control surface. The 
discrepancy is however most likely due to the gap that was introduced between the flap 
and wing in the single grid case to accomodate the shearing grid. If the gap is reduced by 
increasing the spanwise resolution of the wing and control surface, the computational results 
of the refined single grid case (151 x 87 x 34) approach those of the virtual zone. This 
particular example demonstrates the importance of simulating the geometry of the control 
surface/wing configuration accurately. 

The final figure shows the upper surface pressures as well as the instantaneous particle 
traces emanating from the leading edge of the wing (yellow traces), the lower edges of the 
wing at the control surface cut (blue traces), and the upper edges of the control flaps (red 
traces). The full paper will cover the flow physics near the wing/flap junctures in more detail. 

Conclusions 

An unsteady interface algorithm based on the idea of virtual zones has been developed 
for a finite difference code, ENSAERO. The new ‘virtual’ zoning technique simplifies zoning 
complex geometries such as control flaps and makes possible the use of standard multizonal 
codes for configurations difficult to do before. A fast search routine based on a window 
clipping algorithm has also been developed and is sufficiently fast so that new interpolation 
coefficients can be recomputed at every time step. 

Both of the above developments have made practical a complete unsteady Navier-Stokes 
simulation of a forced oscillating control surface on a clipped delta wing in the transonic flow 
regime. The method has been validated against experimental data as well as numerical sim- 
ulation based on a shearing single zone grid. The numerical result confirms that the accurate 
representation of geometry by ‘virtual zones’ is superior to the single zone computations at 
the same grid size. 
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Fig. 1. Planform of clipped delta wing with trailing edge flap. 

Fig. 2. Surface grids and zonal boundaries of wing/flap configuration. 

Fig. 3. Perspective view of trailing edge flap and wing. 

Fig. 4. Perspective view of trailing edge flap and wing with the inboard wing and flap 
virtual zones. 
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and flap virtual zones. 

Fig. 6. Comparison of unsteady pressures between virtual zones and single grids com- 
putations with the experimental data of Ref. 11. 

Fig. 7. Upper surface pressures and instantaneous particle traces emanating from the 
leading edge and the edges of the control surface. 
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Canard- Wing-Body Ramp Motion 

Eugene L. Tu,* Shigeru Obayashiy and Guru P. Guruswamvi 
NASA Ames Research Center, Moffett Field, California 94035-1000 

Introduction 

The use of canards in many advanced aircraft for control and improved aerodynamic 
penormance has been the topic of continued research. In addition to providing rapid pitch 
control, the influence of canards on wing aerodynamics can often result in increased maximum 
lift and decreased trim drag. There are also unique dynamic performance characteristics 
for canard-configured aircraft coupled with the capability of present-day automatic control 
systems. The reduced or even negative static stability of canard configurations can lead to 
improved aircraft agility and maneuverability. 

Several examples of the use of canards ior stability and control are currentlv available. 
The X-31 aircraft uses a long-coupled canard for pitch control 1 while the SAAB JAS 39 
Gripen uses a short- (or close-) coupled canard in maneuvering, cruise and even landing roll- 
out conditions. 2 The close-coupled canard of the X-29 forward-swept aircraft is integrated 
into the active control system and is used to maintain control of this inherently unstable 
aircraft. 3 

In the three examples given above and, indeed, in most canard-configured aircraft, 
the main benefits of canards are realized during maneuver or other dynamic conditions. 
Therefore, the detailed study of canards as primary control surfaces requires the accurate 
prediction of the unsteady aerodynamics of such configurations. For close-coupled canards, 
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the unsteady aerodynamic performance associated with the canard-wing interaction is of 
particular interest. 

At moderate angles of attack, canards or wings with sharp leading edges exhibit flow 
separation at the leading edge due to the adverse pressure gradient on the leeward side. In 
general, the flow structure of highly-swept or delta canard-wing configurations is character- 
ized by a canard downwash, which modifies the wing flowfield, and an interaction between 
the canard and wing vortex systems. The inboard wing flowfield is often dominated by the 
canard downwash and the outboard is affected by the subsequent change in wing leading- 
edge vortex formation and the canard-wing vortex interaction. Further details of the flow 
features of steady canard- wing- body aerodynamics are given in Ref. 4. 

In general, the characteristics of static canard configurations are adequately represented 
by steady-state aerodynamics. At higher angles of attack, some of the conditions which 
may result in unsteady aerodynamics include large regions of separated flow and vortex 
breakdown. However, for a configuration undergoing unsteady motion, the dynamic effects 
can be quite significant. In particular, the downwash of the canard and the interaction 
between the canard and wing vortices can exhibit highly non-linear unsteady aerodynamic 
characteristics. 

The use of canards for improved cruise performance has been supported by both exper- 
imental (e.g. Refs. 5-9) and computational studies (e.g. Refs. 10-12) which investigated the 
steady-state aerodynamics of typical canard configurations. Many of these studies investi- 
gated the effects of canard size, position and deflection angle on the canard- wing aerodynamic 
interaction. A study by Bovden 13 investigated the dynamic stability and response charac- 
teristics of typical canard configurations and showed potential benefits in maneuverability 
and agility with the use of canards. 

Computational fluid dynamics (CFD) has become a valuable tool for understanding 
the complex three-dimensional flow physics of canard configurations. A number of studies 
based on conformal mapping, linear and nonlinear vortex lattice methods, the transonic 
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small perturbation (TSP) equation, and Euler equations have been performed for steady 
canard-wing aerodynamics and are listed in Ref. 4. However, very limited computational 
work has been performed to study the details of the unsteady canard-wing-body flowfield 
using the Navier-Stokes equations, which are required to accurately model viscous effects. 

Previous studies 4 ’ 14 - 15 by the first author have successfully solved the Navier-Stokes 
equations to investigate steady-state canard-wing-bodv aerodynamics, including effects of 
canard deflection and vertical position. Accuracy was demonstrated by favorable compar- 
isons with experimental surface pressure, force and moment data. A grid refinement study 
was also performed to resolve any significant discrepancies in the baseline computational 
comparisons with experiment. 

A more recent Navier-Stokes simulation 16 has been performed to investigate the un- 
steady aerodynamics of a wing-bodv (no canard) configuration undergoing ramp motions. 
Reference 16 demonstrated significant dynamic effects on wing-body aerodynamic loads. In 
the present study, the thin-layer Navier-Stokes equations are solved for the unsteady flow 
about a highly-swept canard-wing-bodv configuration undergoing ramp (pitch-up) motions. 
Emphasis is placed on understanding the complex unsteady flowfield at various pitch rates 
and flow conditions. In addition to surface pressures, forces and moments, a detailed analy- 
sis of the unsteady canard-wing vortex structure, including vortex breakdown, is performed. 
Both grid and time-step refinement studies are conducted to verify adequate spatial and time 
accuracy of the current method. 

Computational Modeling 

Numerical Procedure 


The NASA Ames ENSAERO code is used, to solve the unsteady thin-laver Navier- 
Stokes equations. ENSAER.0 has the capability to simultaneously integrate the Navier- 
Stokes equations coupled with the modal structural equations of motion and has been recently 
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demonstrated for steady, unsteady and aeroelastic applications. 16 19 Since the current study 
is restricted to rigid-body motions, the modal structural equations are noc being solved. 

The current version of ENSAERO includes both a central-difference scheme, 20 which 
is identical to that used in the Transonic Navier-Stokes (TNS) code, 4,21 and a streamwise 
upwind numerical scheme. 18 Both schemes have been found to be accurate for the current 
geometry and flow conditions. In order to compare current results directly with earlier 
canard-wing-body steady-state computations, 4 the central-difference scheme is utilized for 
this study. It is noted that the central-difference scheme in ENSAERO is first-order accurate 
in time and second-order accurate in space. 

To provide turbulence closure, the Baldwin-Lomax algebraic eddy-viscosity model 22 is 
used. Due to the vortex-dominated flow structures of the highly-swept sharp-leading-edge 
canard and wing, a modification to the original Baldwin-Lomax formulation is required. 
For this study, the Degani-Schiff modification, 23 as originally developed for crossflow-type 
separations, is employed. It is noted that with presem-dav CFD technology, higher-order 
eddy- viscosity models could easily be utilized and are readily available within the ENSAERO 
code. However, with the lack of significant non-equilibrium or streamwise separation effects 
anticipated at the moderate angles of attack being investigated, the benefits for the current 
study of such higher-order models do not justify the increased computational costs. Further 
details about the ENSAERO code, algorithm, zonal approach, and general performance will 
be given in the full paper. 

Geometry Modeling and Grid Generation 


The geometry in this study is based on the wind-tunnel model used by Gloss and 
Washburn 5 and is illustrated in Fig. 1. The same geometry was used in the previous steady- 
state numerical study 4 and in the unsteady wing-body (canard-off) study. 16 In the original 
wind-tunnel model, fairings were used to facilitate a vertical-offset canard. These fairings, 
which account for slight asymmetries in the experimental results, are omitted in the current 
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computational modeling. The sting used for wind-tunnel mounting is modeled by extending 
the body, with its appropriate no-slip boundar;/ condition, to the downstream boundary. 

Using the S3D surface geometry and grid generation code, 24 the canard, wing and body 
component surface geometries are modeled from their original analytical definitions. The 
3DGRAPE 2 ° program is then used to generate the single-zone canard-wing-body flowfield 
grid. The overall H-0 topology grid, with the mismatched interface, for a typical high- 
canard case is given in Fig. 2. An expanded view of the flowfield grid near the geometry is 
also shown. 

The resulting grid for the canard-wing-body configuration contains 4,625 points (half- 
body) on the surface, and approximately 470,000 points in the flowfield. This same grid was 
used extensively in previous steady-state computations 4 and was found to provide accurate 
results at moderate angles-of-attack. A grid refinement study which increases the total 
number of flowfield points to almost 2 million will also be presented in the full paper. 

Since the current computations are performed in the transonic regime, the flowfield grid 
is extended upstream and downstream by approximately eight wing root-chord lengths, and 
in the radial direction by six wing-span lengths. 

Results and Discussion 


Preliminary results are presented to validate the current methodology and provide initial 
analyses of the unsteady canard-wing-body flowfield. All steady and unsteady resuits are 
computed at M 0 0 = 0.90 and Reynolds number, based on the mean aerodynamic chord of 
the wing (Reg) of 1.52 million. Comparisons between computed results and experimental 
data 0,26 are made to validate the accurate prediction of the steady flowfield. Convergence 
of the unsteady flowfield is verified using time-step refinement. Previous studies 16,1 ‘ ,2 ‘ have 
demonstrated accurate unsteady flowfield predictions- using ENSAERO on various wing, 
wing-body and wing with control surface geometries. 
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Experimental Comparisons 


Figure 3 shows the effect of the canard on che wing flowfield at 0 ^ 4 °. The stations 
inboard of the canard-tip span line (25% and 45% span) are most significantly influenced by 
the presence of the canard. For the canard-off case, the suction peak on the upper surface 
identifying the presence of a leading-edge vortex is clearly evident. At this static angle 
of attack, the canard-on results show that the formation of the wing leading-edge vortex 
is inhibited at the inner stations. This effect of the canard on the leading-edge vortex is 
directly attributed to the canard downwash. The results presented in Ref. 4 indicate that 
discrepancies between the computational and experimental results in Fig. 3 are resolved 
with grid refinement. 

Effect of Canard on Unsteady Aerodynamic Loads 


The effect of the canard on the unsteady aerodynamic loads associated with pitch- 
up ramp motion are presented in this section. All unsteady canard- wing-body results are 
computed for a from 0° to 15° and various non-dimensional pitch rates as defined by A = 
ac/Uoo- Ramp motions are started from converged steady-state solutions at an initial angle 
of attack (a,) and held at the final angle (a/) for a specified length of time. The pitch axis 
of the ramp motion and the pitching moment results are taken from the model c.g. location 
shown in Fig. 1. 

A typical ramp motion from a; = 0° to a/ = 15° is illustrated in Fig. 4. In order to 
directly compare unsteady and steady-state results, time (t) is given in degrees. Therefore, 
during the ramp motion (0° < t < 15°) a and t are equal (a = t). However, for t > 15°, a 
is held constant at 15° (a = 15° or a/). In this manner, the ramp motion has an impulsive 
start and finish which significantly influence the unsteady results. 

Computed time histories of lift, drag and pitching moments for the configuration, with 
and without canard, undergoing a ramp motion (pitch rate, A = 0.10) are illustrated in Fig. 5. 
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Computational and experimental steady-state results axe also given at the corresponding 
instantaneous angles of attack. 

The lift curves of Fig. 5 show a significant dynamic lift increase for the unsteady canard- 
on case. Early in the motion ft < 4°). the unsteady canard-on lift is slightly lower than the 
canard-off lift due to the effects of the virtual (or apparent) mass of the fluid. Since the 
entire canard is forward of the pitch axis, there is an initial loss of lift on the canard at 
the start of the ramp motion (t = 0°). As the ramp motion continues, the canard-on case 
exhibits increased dynamic lift over both the canard-off and steady-state canard-on cases. 

As expected, the time histories of drag coefficients follow a similar trend. However, by 
replotting the drag results of Fig. 5 in traditional drag polar form, Fig. 6 shows that the 
unsteady canard-on case exhibits improved dynamic lift-to-drag performance. 

The pitching moments given in Fig. 5 also illustrate the significant influence of the fluid 
virtual mass at the start (t = 0°) and finish (t = 15°) of the ramp motion. The virtual 
mass acts to counter the acceleration of the body and, therefore, causes a rapid nose-down 
pitching moment at t = 0° and a nose-up moment at t = 15°. Beyond t = 15°, the lift, drag 
and pitching moment values converge toward the steady-state a = 15° result. 

The effect of pitch rate on dynamic lift for the canard-on case is given in Fig. 7 for A 
= 0.10 and 0.05. Note that since comparisons are made at instantaneous angles of attack 
during the ramp motion (t in deg.), the physical time (t in sec.) between the two pitch rates 
differ by a factor of two. Figure 7 shows that the dynamic lift of the canard configuration is 
increased at the higher pitch rate throughout the ramp motion. Far enough into the ramp 
motion (i > 10°), there appears to be an approximately linear relation between pitch rate 
and dynamic lift for the current configuration. This observation will be explored further in 
' the full paper with a more extensive parameter study. 

Analysis of Unsteady Aerodynamic Loads 


Better insight into the dynamic loads produced by the ramp motion (Fig. 5) can be 



attained by examining the separate component regions of the geometry. The canard region 
consists of the canard and the body forward of the .wing leading-edge root location (fore- 
body). The wing region consists of the wing and the remaining aft-body (not including the 
sting). 

The time history of the canard and wing region lift contributions for the A = 0.05 ramp 
motion are illustrated in Fig. S. Computed steady-state lift coefficients are also given for 
reference. By definition, the total configuration lift is the sum of the canard and wing region 
lift. The effects of the fluid virtual mass is evident at both a,- = 0° and oc/ = 15°. Due to 
the relative locations of the canard and wing to the pitch axis, there is an initial increase in 
wing region lift and a decrease in canard region lift at a,-. These trends are then reversed at 
a/. 

Figure S shows that there is a net loss of dynamic lift for the canard throughout the 
ramp motion. The increased dynamic lift for the total configuration is due to the large 
increase in dynamic lift of the wing. In linear stability theorjr, the lift of the wing or canard 
can be written as 


Cl = Ci a o. -f C^a + C^q (1) 

where CL a ,CL A and Cl^ are the stability derivatives (9Ci/5a)o, (<9Cz,/<3a;)o and (<3C£,/3<?)oi 
respectively; q is the angular velocity of the configuration about the pitch axis. The notation 
()o indicates that partial derivatives are evaluated assuming no disturbance from the other 
terms. Note that for the ramp motion, q = a. 

The loss of canard lift is primarily due to the location of the canard forward of the pitch 
axis which, if described bv linear stabilitv theorv, causes a negative contribution to Cl 
from the C L q q term in Eq. 1. After the ramp motion stops (f > 15°), the lift contribution 
from the wing converges to the steady-state result much more slowly than the canard region 
lift. Previous studies 16 indicate that vortex interaction and breakdown may be a significant 
factor in this phenomenon. 
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In order to confirm the convergence of the unsteady computations. Fig. 9 shows the 
effect of time-step size (A t in deg.) on the unsteady lift curves for the component regions of 
the geometry. From Fig. 9, it is clear that with decreasing time-step size, the time-accurate 
solutions converge quickly. In fact, even when using a larger step size, the lift curves compare 
favorably for t > 0.1°. In till previous computations presented in this study, A t = 0.0025° 
was used. 

A complete analysis of the unsteady canaxd-wing-body flowfield for the configuration 
undergoing ramp motions will be presented in the full paper. This analysis will include 
a detailed study of the canard-wing vortex interaction and the dynamic effects associated 
with vortex breakdown. For example, a side-by-side comparison of upper surface pressures 
(Fig. 10) shows a dramatic increase in canard and wing vortex strengths for the unsteady 
case. In Fig. 10, the unsteady result is for A = 0.10 and gives an instantaneous map of 
surface pressures at a = 8.55° during a pitch-up ramp motion from 0° to 15°. Numerical 
issues including the effects of dissipation, turbulence modeling and grid fineness will also be 
addressed in the full paper. 
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Navier-Stokes Computations on Full-Span Wing-Body Configuration 
with Oscillating Control Surfaces 

Shigeru Obayashi,* Ing-Tsau Chiu* and Guru P. Guruswamy* 

NASA Ames Research Center, Moffett Field, California 

Introduction 

Aircraft are often subject to aeroelastic oscillation, especially in the transonic regime, 
because of flow unsteadiness in the presence of the moving shock waves. In this unsteady 
aerodynamics environment, aircraft rely heavily on active controls for safe and steady flight 
operation. Active control is also needed for the suppression of flutter without adding 
structural weight to achieve stable flight conditions. 

The influence of control surfaces on both aerodynamics and aeroelastic performance of 
a wing is more pronounced in the transonic regime. These influences can be constructively 
used to improve the wing performance through proper actuation of the active control 
surfaces. Active control technology relies on accurate predictions of unsteady aerodynamics 
and aeroelastic performance of a wing. Since the experimental evaluation of the effect of a 
control surface on the wing performance would involve considerable cost and the risk of 
structural damage in a wind tunnel, it is necessary to initiate the investigation through 
theoretical analyses. 

The present investigation is initiated in conjunction with a recently developed code, 
ENSAERO, which is capable of computing aeroelastic responses by simultaneously 
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integrating the Euler/Navier-Stokes equations and the modal structural equations of motion 
using aeroelastically adaptive dynamic grids . 1 * 4 The code has been applied to transonic 
flows from small to moderately large angles of attack for fighter wings undergoing 
unsteady motions. Furthermore, the geometric capability of the code was extended to 
simulate unsteady flows over a rigid wing with an oscillating trailing-edge flap. 

This paper reports unsteady Navier-Stokes simulations of transonic flows over a rigid 
arrow-wing body configuration of supersonic transport-type aircraft with oscillating control 
surfaces. The base grid was generated by using the hyperbolic grid generator . 5 For 
unsteady computations, the grid moves every time step following the deflection of the 
control surface. Computations have been made at moderate angles of attack with and 
without control surface deflections. The flow condition selected is in the transonic regime 
with a moving shock wave, including leading-edge separation. Computed pressures have 
been compared with the wind-tunnel experiment 6 In the full paper, dynamic stability of the 
model will be discussed in detail by using computed results. Comparison of response 
characteristics between symmetric and antisymmetric control surface motions on the right 
and left wings will also be studied. 


Numerical Method 

The nondimensionalized Reynolds-averaged thin-layer Navier-Stokes equations are 
used in this study. The viscosity coefficient is computed as the sum of the laminar and 
turbulent viscosity coefficients where the laminar viscosity is taken from the free stream 
laminar viscosity, assumed to be constant for transonic flows. As an option, Sutherland's 
law can be used to calculate the laminar viscosity. The turbulent viscosity is evaluated by 
the Baldwin-Lomax algebraic eddy-viscosity model . 7 Since the flow field to be considered 
in this paper contains leading-edge separation, it is important to apply a modification to the 
turbulence model originally developed for crossflow-type separation . 8 


Several numerical schemes have been developed to solve the Navier-Stokes equations. 
The present code has two different schemes for the inviscid term: the central-difference and 
streamwise upwind schemes. A second-order central-difference evaluation is applied to the 
viscous term. An implicit method is used for the time integration because it is more suitable 
for expensive unsteady viscous calculations. The complete algorithm can be found in Ref. 
2. Specific code performance information for the current study is given as follows. All 
results were computed on either a CRAY-YMP or CRAY-2 computer at NASA Ames 
Research Center. The performance of the central difference version of ENSAERO is 160 
MFLOPS and 15 jisec per iteration per grid point on a single CRAY-YMP processor. 

Sample Results 

The first sample result is to show the code capability for computing an oscillating 
control surface. Figure 1 shows the planform of the clipped delta wing with the trailing- 
edge control surface. The wing has a leading-edge sweep angle of 50.4 deg and a 6%-thick 
circular-arc airfoil section. At Moo = 0.9 and a = 3 deg, both a leading-edge vortex and a 

shock wave are present on the upper surface of the wing. The present C-H grid contains 
151 x 44 x 34 points. The gap regions are introduced at the both ends of the control surface 
(see Fig. 2). This region is used to shear the grid when the control surface oscillates. 
Although the gap is introduced to simplify the calculations, its effect can be minimized by 
clustering the grid in this region. The dynamic grid around a deflected control surface was 
obtained by shearing every grid line normal to the control surface with the local deflection, 
Ax and Az. Figure 3 shows the unsteady pressures with the control surface oscillating at a 
frequency of 8 Hz and an amplitude of 6.65 deg at Moo = 0.9, a = 3 deg and Re c = 17xl0 6 

based on the root chord. Results are shown as magnitude and phase angle of the upper 
surface pressure responses at three spanwise sections. The magnitude part of the unsteady 
pressures shows significant influence of the control surface oscillation. Overall, the 
computed results show reasonably good agreement with the experiment. 
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The H-H topology grid is used for a wing-body configuration with a control surface. 
The ICEM DDN CAD software system by CDC was used to generate the surface grid. 
Then the volume grid was generated by using HYPGEN code. 5 To treat the control surface 
movement, a small gap is introduced at the end of the control surface. Figure 4 shows the 
geometry of the wind tunnel model. 6 The configuration is a thin, low aspect ratio, highly 
swept wing mounted below the centerline of a slender body. The wing is flat with a 
rounded leading edge. Figure 5 shows the surface grid for the half-span configuration. The 
body is extended to downstream. This half-span grid consists of 110 points in the 
streamwise direction, 116 points in the spanwise direction, and 40 points normal to the 
body surface, in total of 510,400 points. In the following computations, the grid is further 
divided into the upper and lower grids at the wing and the H-topology cut condition is 
provided through a zonal interface. For the full-span configuration, the grid is mirrored to 
the other side and thus the number of the grid points is doubled. It should be noted that the 
exact wing tip definition was not available so the tip thickness was decreased to zero across 
three grid points. 

Figure 6 shows the steady pressures compared with experiment at three spanwise 
sections for the half-span configuration. The flow conditions consists of Moo = 0.85, a = 8 
deg, 8=0 deg and Re z = 9.5 x 10 6 based on the mean aerodynamic chord. No data 
correction was applied to either the computed or measured data. The leading-edge vortex is 
captured by the computed results well. The grid refinement study for the half-span 
configuration will be included in the final paper. Figure 7 illustrates the instantaneous 
antisymmetric position of oscillating control surfaces under computation. Since no 
unsteady measurement is available, computed mean pressures are validated against the 
steady data. Comparison of response characteristics between symmetric and antisymmetric 
control surface motions on the right and left wings will be studied. In the full paper, 
dynamic stability of the model will be discussed in detail by using computed results. 
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Fig. 1 Planform of a clipped delta wing. 
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Fig. 2 Grid for a wing with a control surface. 
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Fig. 5 Surface grid for the half-span configuration. 
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Fig. 7 Instantaneous view of antisymmetric oscillations of control surfaces. 
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